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ABSTRACT 



Fourier transform techniques have been the favored methods in the analysis of 
signals and systems. One major drawback of Fourier methods is the difficulty in 
analyzing transient and/or non-stationary behavior. Recent advances in the field of 
wavelet theory show much promise in alleviating these problems. This thesis considers 
the realizations of the wavelet decomposition and reconstruction algorithms for the 
discrete case. The major discussion will involve both the one and two dimensional 
transforms. We also present a multiple-phase development as a second and possibly a 
preferable method for decomposing signals. 
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I. INTRODUCTION 



The recent advances of wavelet theory have spumed much interest in many 
disciplines: earthquake detection, medical EKG’s, and multiple resolution signal 
processing. We are most interested in the last discipline listed, multiple resolution signal 
processing. Using wavelet theory in this field opens several areas of interest. In the 
one-dimensional case, analysis of short-time duration signals such as radar pulses, can 
possibly lead to better discrimination of similar but distinct radar sources. For the two- 
dimensional case, statistical pattern recognition using wavelets can aid in identifying 
skewed or rotated images. Lastly, for both dimensions, wavelet theory can be used in 
data compression, which has major applications particularly in transmitting data on a 
bandwidth-constrained channel. Traditionally, most means to date use Fourier methods 
to attack the problems in the areas of interest discussed previously. However, for 
transient and/or non-stationary phenomena, Fourier techniques give inadequate results in 
the transform domain, even when we modify them to include a dimension of time. Thus, 
wavelets can serve as an alternative to the conventional Fourier transform techniques for 
analyzing such transient and/or non-stationary behavior. 

This thesis develops one and two-dimensional invertible discrete wavelet transforms 
(DWT). 386MATLAB and PROMATLAB Version 3.5 are the software tools used for 
386 or better PC’s, and for the UNIX-based SUN workstations. We did not use Fortran 
and C programming languages to maximize the portability of the routines, but such codes 
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could be developed very easily from this thesis. These routines can be used by other 
researchers for more in-depth analysis of signals and systems. 

Chapter II will cover the basic principles of the wavelet theory, particularly for the 
discrete dyadic case. Although many volumes are dedicated to the field, this cursory 
outline of the theory will provide the basis for notation used throughout the thesis. 
Chapter III will entail some important concepts such as causality and spillover effects that 
must be accounted for in the DWT. We expand these ideas further into the two- 
dimensional case in Chapter IV. Chapter V proposes another method for decomposing 
the data at the expense of physical memory, the so-called multiple-phase DWT 
(MP/DWT). The MP/DWT possibly approximates the continuous WT, with the property 
of time invariance. Chapter VI gives a narrative description of the MATLAB algorithms 
developed. Finally, we present our conclusions in Chapter VII. 
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n. BASIC WAVELET THEORY 



Since wavelets are a relatively new and less established field to most readers, one 
main question asked is "what are they?" We introduce the topic in here by an example. 

Consider a function f(t) that has discrete levels as shown in Figure 1. We 
decompose f(x) into a lower resolution level. In other words, we desire to smooth (or 
average) out f(x), or to low pass filter the function to a coarser level. Figure 2 depicts 
the resulting decomposition function „f(x), where a is a scaling factor. For the dyadic 
case, we set a to two. If we compare the original function f(x) to the coarser function 
„f(x), we note that the detail is „d(x)=f(x)-^f(x) as shown in Figure 3. If we had the 
detail and the decomposed signal, we can theoretically reconstruct f(x) without any loss 
of information. The decomposition (and the reconstruction) process can go to any 
arbitrary resolution level m. Usually, we sample the signal at a rate higher than or equal 
to the Nyquist rate, and set the resulting sampled data to be at the highest resolution level 
(say level m=0). All other resolution levels would be in the negative direction. This 
process could be repeated for the next lower resolution level by operating on the current 
level in an iterated manner. 
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Figure 1. Example Function f(x) 
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Figure 2. Blurred „f(x) Figure 3. Detail „d(x) 
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The primary goal is to find a set of orthonormal basis functions that will do the 
smoothing of the original level and retain the details for reconstruction in a more 
practical manner. The low pass filtering process uses a function called the scaling 
function 0(x) while the detail basis function uses the wavelet basis function \p(x), both 
of which are also orthogonal to each other. 

The described method is a coarse level in the understanding of the wavelet theory. 
We must now go to a higher resolution level of knowledge to better comprehend the 
following chapters. 

A. THE SCALING FUNCTION <j>{x) 

The region L^(R) depicts a vector space where a function f(x) resides and satisfies 
the condition of having finite energy, or the inner product of f(x) with its complex 
conjugate is finite. Consider a family of embedded closed subspaces, V,„, such that the 
infinitely largest subspace is L^(R) while the infinitely smallest subspace {0}. 

■fop 

/ f(x)f(x)dx = {/(x),/*(x) > < 00 (1) 

— OP 

(lower resolution ... C V., C Vq C V+, C ... higher resolution) 
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Additionally, for each vector space in let W„ be an orthogonal complementary 
space to in such that V„+,=V„0Wn,. We show the spaces graphically in 
Figure 4. 




Note that if given V„ and to we can obtain the vector space as shown 

in Figure 5. 



V 




m m+1 




V, 



m+k+l 



Figure 5. Another View 
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A function f(x) will lie in a vector space V„, if and only if it can be written as: 



m = E ® 

n € Z 

where we adjust the width of <^(x) for unit grid spacing at resolution level m=0. The 
inner product of 4 > with an n-shifted version of itself must equal zero 
( <</)(x),</)(x-n)> = 0 ). 

In order to normalize the scaling function for the corresponding resolution level in 
the L^(R) space, multiply the scaling function by a factor 2“^^ while also dilating in the 
spatial domain. The following two equations simplify the notation for the rest of the 
derivations. 



= 2hc(2"x-n) = 4>,(x-2-n) 

Let P,n be the orthogonal projection of a function in L^(R) onto the vector space 
V„,. Note that Equation 2 is a harmonic series representation of the function projected 
onto the V„ space, with Cn„, denoting the weighting coefficients for the basis function. 
This is similar to the evaluation of the fourier coefficients for the fourier series expansion 
with its basis functions. 



7 



Thus Equation 5 can be viewed as a convolution of the function with the mirrored 
scaling function sampled every 2” seconds. 

= /W ♦ l2-„ 

The recursion property for the scaling function only depends on two adjacent resolution 
levels m and m+1. Since V„ C V„+j, we have: 

^ ^ (7) 

k 

We define the filter coefficients in terms of the inner product of Equation 7 as 

h(K) . f^,(x) . 2 '^ ( 4i„ , 4>„ ) ® 

We are now able to bring the previous equations into a form known in the wavelet 
field as the fundamental dilation equation. This is a two scale difference equation 
relative to the two adjacent levels. 

4>(x) = 52 2 h(k) ^{2x-k) (9) 

k 

A finite set of h(n) coefficients gives very desirable traits for selecting wavelets 
with compact support. Major research in the field has been in the study and 
determination of the family of scaling functions that are of compact support and 
orthonormal. With the scaling function being orthonormal, we define the following 
properties which are relative easy to prove [Ref. 1: pp. 689-691]: 
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1 <Kx) dx = 1 


(10) 


E m = 1 

k 


(11) 



E = j (12) 

k ^ 



2 Hk) Kk-ln) = 6(«) 
k 



(13) 



Going into the spectral domain, the low pass properties of the scaling function are 
obvious from the next four equations. 







where H(f) = h(k)e-^'^^^ 

it 



(14) 



^(f)= n 

p=i 



H 



if] 

W) 



(15) 



Y, Hk) = 1 - HiO) = 1 (16) 

it 

I H(f) 1^ * I /f(/4) p = 1 (17) 
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B. THE WAVELET BASIS FUNCTION ^(x) 

In a manner analogous to the development of the scaling function, the wavelet basis 
function \p(x) determines the W„ vector space. In a similar fashion to Equation 2, we 
define the detail function as follows; 

m 

= E 2’' ^(2’x-n) (18) 

neZ 

The wavelet basis function is of the form of a constant-Q band pass filter function, 
increasing an octave in bandwidth as the frequency increases an octave. Since is 
orthogonal to V„, the inner product of each of the spaces’ basis functions should equal 
zero ( < \p(x) , 4>(x) > = 0 ). This relationship is so intertwined that by knowing a 
valid set of a set of can be derived for a particular resolution level. 

Consequently, knowing the h(n) coefficients will determine the appropriate coefficients 
for the wavelet basis functions given as g(n). 

We define the wavelet basis in a similar fashion as the dilation property in Equation 
9: 

^{x) = 2 g{k) (p{2x-k) ( 19 ) 

k 

The {^m(x-n)} is an orthonormal set spanning W„. The g(n) coefficients share many of 
the properties of the h(n) sequence with the following exceptions: 
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i'{x)dx =J^^(^)=0 

k 



( 20 ) 



m 




( 21 ) 



Mallat [Ref, 1: p. 681] and Daubechies [Ref. 2: p. 944] chose the g(n)’s in terms 
of the h(n)’s as follows: 



gin) = (-1)*-" h{\-n) 



( 22 ) 



C. THE DISCRETE WAVELET TRANSFORM 

Note that Equation 6 shows the lower resolution, or the approximation, signal’s 
weighting coefficients as a convolution of the time reversal of the h coefficients with the 
signal f(x) sampled every 2“ seconds. Using the recursive properties of the scaling 
function, it is easily shown that the approximation coefficients of the lower resolution 
levels also can be determined from approximation coefficients of the higher level. In 
particular: 
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^mn ^ > ^mn ^ 

= 1 <i>m.ldX 

“E ^ ^m’^/wn+l^ 

k 

=E 2’^ c„., „ 

k 

=2^5^ *(t) c„.,(2n*t) 

Jt 

=2’5:^(2«-;X.,0) 



where h(n) = h(-n). 



The detail weighting coefficients can similarly be determined as 



- 2 " T,lQn-j)cJ!) 



Recall that the d^’s are the detail coefficients of the orthogonal projection of f(x) onto 



the W„ vector space. The entire collection of {d^}, or {c.^n} U {dnm • > -M} is 



called the discrete wavelet transform. 



1. Decomposition 



Decomposing to any lower level requires that the L^(R) function be initially 



convolved with <^o(-x), and sampled at a unity interval to obtain the coefficients, Co(n), 



as shown in Equation 2, and graphically depicted next: 
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c 



On 



Figure 6. Determining m=0 Coefficients 

Once the Con coefficients are known, we use Equations 23 and 24 in an 
iterated manner to obtain the set of coefficients for any arbitrary level. For the dyadic 
case, this process is a convolution, a two-times decimation of the result, followed by a 
2’^^ scaling factor. 




m-l,n 



Figure 7. Decomposition Algorithm 



13 



2. Reconstruction 



Recall that V„+i=V„© W^. The c^’s reside in the space and the d^’s 
lie in the space. Thus, a projection of a function onto the next higher level should 
equal the sum of the projections onto V„, and W„. Let and Qm be the projection 
operators onto the vector space. 

_ (25) 

n 

QJfi = 

n 

Matching the terms in Equation set 25 and doing several transformations of 
variables, we get the following: 

‘^..1 E \cJl)Kn-2l) (26) 

I 

Equation 26 shows that the reconstruction of the approximation coefficients 
of one level higher can be realized by putting zeros between the lower resolution 
coefficients, convolving with the respective H or G filter, adding and then scaling the 
result by 2^'^. We show a block diagram of this basic algorithm in Figure 8. 
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c 

ni-l,n 



d 

m-l,n 



Figure 8. Reconstruction Algorithm 

For the rest of the discussion, we will concentrate primarily on the 
approximation and detail coefficients. It is the properties of these coefficients that are 
of importance in relation to the DWT process. 

The basic procedure would be to decompose the data array, starting from the 
highest resolution level. Next calculate the coefficients for the next lower level, while 
retaining the details, d„u„ and using the resulting approximation coefficients to calculate 
the next lower level. Repeat this procedure until reaching the lowest desired level where 
now both the approximation and detail coefficients are retained. 
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m. ONE DIMENSIONAL DWT DEVELOPMENT 



A. DECOMPOSITION 
1. Causality 

Recalling that the decomposition algorithm of Figure 7, a finite data vector 
of the sampled signal, Con, can be decomposed to its lower resolution coefficients, c.,„ and 
d.in. This process requires that the data array be convolved with the respective filters. 
These filters are the time-reversed version of the H and G vector array. At the output 
of each convolution, there is a two times decimation. The scalar multiplication 
normalizes energy the in the L^(R) domain. 

Consider a finite array Co„ of length L. Let Con to be a non-zero set only if 
the indices, n, satisfy 0 < n < L-1. Also set the number of h coefficients equal to an 
even number N. By using Equation 23, the convolution process can be depicted as 
shifting the h’s by a factor of two each time. 



o 

o 

o 


o 

o 






o 

O 


o 

o 

ih 




Co., 


O 

o 


c 

'^0,10 



Shift by two 







Shift by two 



Figure 9. Shift Factor of 2 for N=4 
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In an effort to make the transform a causal system, we desire to make all the 
non-zero data coefficients have indices greater or equal to zero since negative indices 
would equate to negative time. This reconstruction gives two possible sets of indices for 
the h filter coefficients. 

• {h(-N-t-l),...,h(-2),h(-l),h(0)} 

• {h(-N4-2),...,h(-l),h(0),h(l)} 

We allow the negative indices for the filter since we already know the filter 
coefficients (a-priori data). These two sets for the filter coefficients equate to the two 
possible phases that can be used for the decomposition. For notation purposes, we will 
set the largest value of the indices to equate to the type of phase decomposition, phase-0 
or phase- 1. 

By using Mallat and Daubechies selection of the g coefficients in Equation 
22, there is a problem that occurs. While the resulting convolution with the h mirror 
filter gives the values for indices that are greater or equal to zero, the output of the detail 
coefficients give values for negative indices. In other words, given we define the h’s 
with phase-0 indices, the resulting g coefficient indices are all greater or equal to one, 
as shown in Equation 23. 
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In order to get both the c„„’s and d^n’s to fall into the right hand side, we 
must find another relation for g in terms of h, that satisfies the properties in Equations 
18-21. The following equation will serve the purpose: 

g(n) = (-l)'^*^/z(-iV+3-n) 

2. Phase Selection 

The decomposition can proceed with the choice of phase to select. The 
selection of the desired decomposition phase and size of the input vector, L, determines 
the total number of coefficients at the lower resolution level. 

phase-0: |{c^,J| = 

(28) 

phase-l: 



N+L-1 



2 

N+L-1 



3. Zero-Padding of the Data Array 

The basic algorithm uses the dot product between the h (and g) coefficient 
vector with a set of data coefficients, Cn„,’s, that are obtained by a sliding window on the 
input coefficients, c„,+i „ vector. The sliding window shifts to the right on a generated 
work vector to simplify the programming. This work vector is a zero-padded version 
of the input coefficient vector. Depending on the length of the input data array, the 
selection of the phase, and the size of the filters, the required zero-padding can be 
determined for the beginning and the end of the sequence. 
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If we select a phase-0 decomposition, we will generate the work array by 
initially padding one zero on the left. Additionally, if the size of the data array is even, 
we append another zero to the end of the sequence. For the phase- 1 case, we pad the 
work array with one zero at the end of the sequence only if the original array is odd. 
After determining all these conditions, we finally pad the work array with N-2 zeros both 
in the beginning and the end of the intermediate sequence. 



Initially pad with one zero if the 
Phase-0 : the number of coefficients in the 

array is EVEN. 




Figure 10. Definition of the Work Array 

In esoteric terms, a phase- 1 decomposition for each level would be more 
efficient. Practically, this savings in memory is very insignificant, since the number of 
practical resolution levels is approximately flogzC I Con | )1 • Additionally, this reduction 
in the size of the workspace is balanced by a larger workspace during the reconstruction 
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process. Most importantly, desiring only to decompose the data with only the phase- 1 
case for each level may miss some properties of a signal, which will be further discussed 
in Chapter V. 

4. Energy Determination 

One quick method in assessing if the decomposition routine is working 
without doing the reconstruction is to check the energy of the signal in each resolution 
level. Recall that Equations 2 and 18 are harmonic series representations. By using 
Poisson’s Summation Formula (PSF), the energy in the function in the spatial domain is 
equal to the total energy stored in the weighting coefficients of the transform domain. 
This property holds true when the basis functions are orthonormal in the L^(R) space. 
Thus, the energy of a higher resolution level’s approximation coefficients is equivalent 
to the sum of the energies of the approximation and detail coefficients one resolution 
level down. 



E = I 

n I 

By normalizing the energy to the input resolution level, m=0, a quick check 
of all the levels will add confidence to the program accuracy. In addition, this energy 
check shows the numeric truncation errors inherent on the processing unit used. The 
display of the energy will provide further insight on possible compression schemes for 
classes of signals, since the energy may be only concentrated at certain resolution levels. 



Cm/ +4/ 



(29) 
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Now as an example, consider a one-dimensional transient signal’, with a 
phase- 1 decomposition for all the levels and using a Daubechies 2-tap filter (N=4) 
[Ref. 2: p. 980]. The following diagrams show the original signal and the energy plot. 



X 1 o — •= 




Figure 11. Transient Signal 



Professor Michael A. Morgan, Chairman of the Department of Electrical and 
Computer Engineering at the Naval Postgraduate School, provided the transient data 
included in the MATLAB routines called "et90.m ." This signal is the back-scattered 
transient electromagnetic field from a 10cm long thin-wire illuminated by a double- 
Gaussian impulse, computed using a time-domain integral equation. The other transient 
signal provided by Professor Morgan is "et60.m". The numbers 90 and 60 correspond 
to the angle of incidence that the double-Gaussian impulse impinges onto the wire’s axis. 
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FR«soluitior» 




Figure 12 . Energy in each Level 



The energy check going up from the lowest resolution is depicted next. Note 



that the Daubechies h filter coefficients are only significant to places, thus they 



cause insignificant numerical errors in the reconstruction. 



Verification of the energy in each resolution level 



Level Energy in 
level m-1 
c+d 



Energy in C Difference 

level m 



-7 0, 000024166621151 

-6 0,000258409042490 

-5 0.001707402900115 

-4 0,015878270678385 

-3 0.197762381708774 

-2 0,474398688169256 

-1 0,936969604304463 

0 1,000000000000806 



0,000024166621151 

0.000258409042489 

0,001707402900114 

0,015878270678372 

0,197762381708634 

0,474398688168927 

0,936969604304102 

1,000000000000000 



0,000000000000000 

0.000000000000000 

0,000000000000001 

0,000000000000013 

0,000000000000140 

0.000000000000329 

0.000000000000361 

0.000000000000806 
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5 . Approximation and Detail Time/Scale Diagrams 

One interesting result of the DWT for the one dimensional case is the display 
of the energy on the time/scale plane. The scale can easily be related to frequency by 
noting that the highest resolution is the sampling frequency at m=0. As the scale 
decreases one level, we halve the center frequency for this scale. 

Figure 7 shows the decimation of the coefficients by a factor of two. When 
we obtain the output of the decomposition process, the coefficient’s are in a packed 
format, meaning that the appropriate time spacing is not preserved. Thus, 2”-l zeros 
must be padded between each coefficient for a particular resolution level m. 

The following time/frequency diagrams show the transient signal, mentioned 
in the previous section, with the proper spacing of the coefficients’ energy in the 
transformed domain. 

Zoorr-iecd Energy Disploy of tHe Ap p ro x 1 nn o 1 1 o n Rhose: 11111111 




Figure 13. Mesh Plot for the Energy of the Approximation 
Coefficients to Eight Resolution Levels 
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2oorm®cJ Contour Diaploy of tHa Ap p ro x J m o t T o r* 



F=> H a a a : 



1111111 




Figure 14. Contour Display for the Energy of the 
Approximation Coefficients to Eight Resolution 
Levels 



Zoomed EnerQy Display of tine Detoll RKnose: 11111111 




Mesh Plot for the Energy of the Detail 
Coefficients to Eight Resolution Levels 
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Figure 15 




I 



Zoo m ecd Oontour Disploy of D«toM RJnose: 11111111 




Figure 16. Contour Display for the Energy of the Detail 
Coefficients to Eight Resolution Levels 



B. RECONSTRUCTION 

Figure 8 and Equation 25 show the realization of the recomposition routine by first 
inserting a zero after each coefficient. Next we convolve the array with their respective 
filter, sum, and finally normalize by a factor of 2*'^. This process is essentially an 
interpolation. 
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Since we have defined the indexing of the h and g coefficients during the 
decomposition phase, the recomposition process is fixed to one of two methods 
depending on the decomposition phase used to get to that resolution level. The same 
dyadic convolutional algorithm can be used for the recomposition after initially massaging 
the input data arrays, and shifting one unit to the right instead of two. Initially, the h’s 
and the g’s are reversed in order. 

Now we generate the workspace vectors for the respective coefficients. Append a 
zero AFTER each of the input coefficient vectors. Then add two more zeros at the end 
of the modified work vectors. Finally, if the phase-1 decomposition process was utilized 
to obtain the coefficients, append another zero at the beginning of the modified sequence. 
The diagram below summarizes the work vector definition prior to running the 
convolution algorithm. 



N-2 



p 




0 




0 




0 




0 




0 




0 



0 




Pad one extra zero for the phase-1 case ONLY 
Pad one zero AFTER every "packed" coefficient 
Individual coefficients from the "packed" format 



Figure 17. One-dimensional Reconstruction Work Array 
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All the information needed for the reconstruction is the selected phase 
decomposition to get the c’s and d’s for that resolution level. One factor to note is that 
with Equations 26, there is a possibility that an extra coefficient be generated in the 
reconstruction process. This extra coefficient will theoretically equal zero, but 
numerically it may not and can be cascaded into a notable error if not accounted for. So 
since we know the length of the input vector at resolution level m=0, the size of the data 
array will be known for all levels since we retain all the detail coefficients. Any extra 
zero-valued coefficient generated can be identified as such and ignored in going up to 
the next resolution level. 

Errors in the regeneration will be totally manifested as numeric errors of the 
computer processors. The following diagram shows the reconstructed transient signal 
and absolute error with the original. 




F? « a o I 1:» o n O 

Figure 18 . Original Transient Signal and Reconstructed 
Versions 
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Magnitude of the Absolute Error of 
Reconstruction. Notice the maximum value. 
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rV. TWO DIMENSIONAL DWT DEVELOPMENT 



A. INTRODUCTION 

Arguments for the DWT can be generated for any higher dimension. Of particular 
concern is the two-dimensional case for image processing. The vector subspaces now 
will reside in the L^(R^) vice L^(R) space. A two-dimensional scaling function <^(x,y) 
exists whose scaled dilations and translations for a particular resolution level form an 
orthonormal basis for the vector subspace V„. By making the two-dimensional scaling 
function separable, <^(x,y)=<^(x)<^(y), each vector space can be decomposed as a tensor 
product of two identical sub spaces in L^(R). 

K = vl® Vl (30) 

m fn m 



The tensor product can be viewed as the in-phase and quadrature-phase components in 
V„. The properties depicted in Figures 4 and 5 still apply for the two-dimensional case. 

V = , <8> 

^m+1 ^m+1 ^ ^m+1 



- K® K)‘^ K® I'l) 

= (kI® vl)e (p1® wl)e (K® il)® K® K) 



(31) 
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Notice that defines a set of orthonormal basis functions in and in W„‘. 
Thus, Equation 31 implies that four sets of equations form the orthonormal basis of V„+, 
we summarize in the following table. 



Basis Vector Space Basis Frequency 

Function Coefficient Characteristics 


0mn(y) 


0 yj 


c„(k,l) 


Vert Low Pass 
Horiz Low Pass 


^(^)nm ^mn(y) 




'd„(k,I) 


Vert Low Pass 
Horiz Band Pass 




0 


^d„,(k,l) 


Vert Band Pass 
Horiz Low Pass 


’/'(X)nm ’/'mn(y) 


w„,‘ 0 




Vert Band Pass 
Horiz Band Pass 



Let P„, and denote the projection of a L^(R^) function f(x,y), onto 

vector subspaces <2) V„„ V„ 0 W„, W„, 0 and 0 Wn, respectively. Then in 
a fashion analogous to Equation 25, the vertical and horizontal low pass operation of 
f(x,y) onto resolution level m is as follows: 







(32) 



Notice that the weighting coefficients, Cn,, are now indexed in two-dimensions. 
Similarly, the and ^d„ coefficients can be determined which correspond to the 

low horizontal and high vertical, high horizontal and low vertical, and high horizontal 
and high vertical frequency details. 
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Since the three two-dimensional wavelets and the one low pass basis functions are 
separable, the two-dimensional case is very easily realized. The weighting coefficients 
can be determined by first decomposing on the rows, then the columns of the array, or 



vice versa. The general algorithm is depicted in Equations 33-36 and Figure 20. 



k I 

(I J) - 2^'£c„^k,l)h,^l-y)g^(.k-2i) (34) 

k / 

= 25]5;c„(*:,0«»(/-2/)A,(i-2i) (35) 

k I 

‘ 2£YlcJkJ)s,(l-V)g,(l‘-20 (36) 

k I 




Figure 20. 2-D Decomposition Routine 
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B. DECOMPOSITION 



We will use the convention in image processing that the upper left comer is the 
starting and the lower right is the ending data point. This way, we extend the idea of 
causality into two dimensions, where the spillover effects will fall to the right and below 
the input array. 

1. Decomposition Mask 

By looking at the ‘d„.i coefficients’ indices between the left and right side of 
the equation, generate a two-dimensional mask by forming the outer product of the 
vertical and horizontal filter coefficients. 

- W • [AJ (37) 

where [h] = [h(-N-f2) ... h(0) h(l)]. Similarly the masks for GhH^, and G^Gv also 
can be easily determined. 

We require four masks to decompose successfully an image as suggested in 
Figure 20 to generate the coefficients c„,.i, 'dn,.i, ^d„,.i, and H^Gv, GhH„, and 

GhGv are the corresponding filter masks as shown next: 
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h(-2)h(-2) 

h V 


h(-l)h(-2) 

h V 


h(0)h(-2) 

h V 


h(l)h(-2) 

h V 


h(-2)h(-l) 

h V 


h(-l)h(-l) 

h V 


h(0)h(-l) 

h V 


h(l)h(-l) 

h V 


h(-2)h(0) 

h V 


h(-l)h(0) 

h V 


h(0)h(0) 

h V 


h(l)h(0) 

h V 


h(-2)h(l) 

h V 


h(-l)h(l) 

h V 


h(0)h(l) 

h V 


h(l)h(l) 

h V 



Figure 21. Mask 



h(-2)g(-2) 

h V 


h(-l)g(-2) 

h V 


h^0)^(-2) 


Wl)g(-2) 

h V 


h(-2)g(-l) 

h V 


h(-l)g(-l) 

h V 


h(0)g(-l) 
n V 


h(l)g(-l) 

h V 


h(-2)g(0) 

h V 


h(-l)g(0) 

h V 


h(0)g(0) 

h V 


h(l)g(0) 

h V 


h(-2)g(l) 

h V 


h(-l)g(l) 

h V 


h(0)g(l) 

h V 


h(l)g(l) 

h V 



Mask 



Figure 22 
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g(-2)h(-2) 

h V 


g(-i)h(-2) 

h V 


g(o)h(-2) 

h V 


g^l)h^-2) 


g(-2)h(-l) 

h V 


g(^-l)hj-l) 


g(0)h(-l) 

h V 


g(l)h(-l) 

h V 


g(-2)h(0) 

h V 


g^-l)h(0) 


g(0)h(0) 

h V 


g(l)h(0) 

h V 


g (0) h (1) 

h V 


g(-l)h(l) 
h V 


g|^0)hjl) 


g(i)h(l) 

h V 



Figure 23. Mask 



g(-2)g(-2) 
n V 


g^(-i)^(-2) 


g(0)g(-2) 

h V 


gii)g(-2) 
n V 


g(-2)g(-l) 

h V 


g(-i)g(-i) 

h V 


g ( 0 ) g ( - 1 ) 
n V 


g (1) g ( -1) 

h V 


g^(- 2 )g 40 ) 


g(-l)g(0) 

h V 


g (0) g (0) 

h V 


g (1) g (0) 

h V 


g(-2)g(0) 

h V 


g(-l)g(l) 

h V 


g (0) g (1) 

h V 


g(l)g(l) 

h V 



Figure 24. Mask GhGv 
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These masks generate the corresponding coefficients by shifting by a factor 
of two to the right and downward through a workspace array. With the one-dimensional 
case, there were only two possible phase decompositions. Now four possible phases 
exist. Obviously, the workspace array is different for each phase decomposition. 




Phase-00 




Phase- 10 




Phase- 10 




Phase- 11 



Figure 25. The Four Possible 2-D Decomposition Phases 
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2. Zero-padding of the Data Array 



As stated from the previous section, the work space arrays are different from 
each other depending on one of tlie four decomposition phases selected, the size of the 
mask, and the size of the input array. These considerations must be noted. 

The generation of the work array for the phase-00 case first consists of 
determining if the input array rows and/or columns are odd or even. If either of them 
are even, append a corresponding extra zeros row or column to the bottom or right of 
the input array, respectively. Given the number of rows and columns for the mask as 
Nv and Nh, add Nv-1 zeros rows above and Ny-2 zeros rows below the modified input 
array. Append Nh-1 zeros rows to the left and Ni,-2 zeros rows to the right of the 
modified input array. 



N 



h 



-2 



Nv- 




If # columns is even 
add one more column 
of zeros 



Figure 26. Phase-00 Work Array 
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For the phase-10 case, if the original array’s rows are even and/or if the 
columns are odd, add one respective zeros row/ column to the corresponding bottom/right 
of the modified array. Next add N^-l zeros rows above and Ny-2 zeros rows below the 
work array. Lastly, add Nh-2 zeros columns to the right and left of the array. 



N 



h 



-2 
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o 


o 


o 


o 


a 




o 


o 


o 


o 


o 



} N„- 



N ^-2 



If # rows is even 
add one more row 
of zeros 



If # columns is odd 
add one more column 
of zeros 



2 



Figure 27. Phase-10 Work Array 



Alternatively the phase-01 case, if the original array’s rows are odd and/or 
if the columns are even, add one respective zeros row/column to the corresponding 
bottom/right of the modified array. Next add N^-2 zeros rows above and below the work 
array. Lastly, add N,,-! zeros columns to the right and Nh-2 zero columns to the left of 
the array. 
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N 



h 



-2 



Nv- 




If # columns is even 
■ add one more column 
of zeros 



Figure 28. Phase-01 Work Array 



Lastly for the phase-11 case, if either of the input array’s rows or columns 
are odd, then include a corresponding extra zeros row or column to the bottom or right 
of the input array, respectively. Given the number of rows and columns for the mask 
as Ny and Nj,, add N^-2 zeros rows above and below the modified input array. Add Nh-2 
zeros rows to the left and to the right of the modified input array. 
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N 



h 



-2 




Figure 29. Phase-11 Work Array 



3. Energy Determination 

Similar to the one-dimensional case, adding up the squares of the four sets 
of coefficients at resolution level m will equal the energy in the coefficients for the next 
higher level, m-t- 1. This is a natural two-dimensional extension of Poisson’s Summation 
Formula, discussed previously in Chapter III. A. 5, as shown as follows; 

EE K..P = EE( * W • W ] 

j k j k 
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As an example, consider an image of a centered square, in Figure 30, 
decomposed down to three resolution levels with a phase-00 HAAR wavelet. The energy 
plot in Figure 31 shows the distribution of the c’s and the four d’s for each level. 



Square Image 




3D display of the Image 




Figure 30. Square Image 
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Figure 31. Energy Distribution 
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Except for numerical computations of the computer processor, all the energy 
from the lower d coefficients plus the energy of the lowest c coefficients add up to the 
energy of the original coefficients, Cq. The following energy check results for the two- 
dimensional case, normalized to Cq. Notice that there are no errors since we used Haar 
h filter coefficients and that the two-dimensional decomposition uses a normalization 
factor of two vice the square root of two. 



Verification of the energy in each resolution level 



Level Energy in 

m level m-1 

c+dl+d2+d3 



Energy in C Difference 

level m 



-2 0.440942796610169 0.440942796610169 0.000000000000000 

-1 0.751059322033898 0,751059322033898 0.000000000000000 

0 1.000000000000000 1,000000000000000 0.000000000000000 



C. RECONSTRUCTION 

Since the scaling and wavelet basis functions are separable, the reconstruction 
algorithm closely follows the same process as in the one-dimensional case. First each 
row/column is zero padded appropriately and reconstructed with the one-dimensional 
technique. Then process each resulting column/row in a similar fashion but with a one 
unit shift to the right and then downward, as in a raster scan. We show this process in 
Figure 32. 
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Figure 32. 2-D Reconstruction 

By taking another look at Figure 32, a more efficient method can be realized. By 
padding each coefficient in the array with one zero in both the horizontal and vertical 
directions, we can then slide a reconstruction mask through the padded array. 

1. Reconstruction Mask 

We generate the reconstruction masks in a similar fashion as in the 
decomposition case, but with the important note that the coefficients are in the reverse 
order to accommodate the reconstruction convolution depicted in Figure 32 and the 
following equation in matrix notation 

GA ■ ■ Kj 

where h = [h(l) h(0) ... h(-N+2)]. 
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Another simpler way to generate these masks is to reverse the order of the 



decomposition mask in both the horizontal and vertical directions. We show the 
comparison of the two masks in Figures 33 and 34 for the case only, all other 
reconstruction masks can be similarly related. 



g ( - 2 ) ^( - 2 ) 
n. ^ 


g (-l)h(-2) 
Ji ■V' 


g^0)h^-2) 


g ^1) h^-2) 


g(-2)h(-l) 




g (0) ( - 1) 

K V 


g(l)hL(-l) 

h V 


cr(-2)hL(0) 




g(o)hL(o) 


g(i)h(0) 

tx V 


g ^0 ) 1 ) 




g^OhJl) 


g (1) hid) 



Figure 33. Decomposition Mask Gi,H„ 



g(l)h(l) 

h V 


g(0)hi(i) 


g(-l)h(l) 

hx V 


g 2 ) l^( 1 ) 


g (1) h(0) 

V 


g ( 0 ) hi ( O ) 

ll V- 


g(-l)h(0) 

lx V 


g(-2)h(0) 

lx V 


g (1) h( - 1) 


g (0) h ( - 1) 

hi V 


g(-l)h(-l) 
ll V 


g(-2)h(-l) 
lx V 


g^l)hd2) 


g(0)h(-2) 

n V 


g 1 ) hj - 2 ) 


g ( - 2) h ( - 2 ) 

lx V 



Figure 34. Reconstruction Mask HyG^ 
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2. Zero-padding of the Coefficient Arrays 

In order to use the reconstruction masks, we generate the work array for the 
lower level coefficients to accommodate for the decomposed coefficient array and mask 
size, and the decomposition phase that was selected. 

The phase-00 reconstruction work array is the smallest of the four work 
spaces. The coefficient array is padded initially with one zero to the right and bottom 
of each of the individual coefficients. This makes the modified array have an even 
number of rows and columns. Then append Nt-2 columns of zeros to the right of the 
array followed by Nv-2 rows of zeros on the bottom. 



N -2 rows of zeros 
h 





® i 






I ^ 




0 ( 


) 0 


0 


0 




® i 


i ® 




I 


0 


0 ( 


D 0 


0 


0 




1 ^ i 


i ® 






ol 


0 ( 


) 0 


0 


0 



0 0 0 0 0 0 
0 0 0 0 0 0 



0 0 
0 0 
0 0 
0 0 
0 0 

0 0 
0 0 
0 0 



N -2 columns of zeros 

V 




- Initial Array Coefficients in a packed format 



Padded zero put to the right, bottom or the bottom right of the 
individual packed coefficients 



Figure 35. Phase-00 Reconstruction Work Array 
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For the phase- 10 case, add one extra column of zeros to the left of the phase- 



00 work array. Alternatively for the phase-01 case, append one extra row of zeros to 
the top of the phase-00 work space. Lastly for the phase- 11 situation, insert both one 
row and one column of zeros to the top and left of the phase-00 case. 



Phase-00 
Reconstruction 
Work Array 

from the previous 



Add one column of zeros for 
either : 

Phase-10 
Phase- 1 1 



Figure 36. Work Space for Phase-10,10,11 Cases 
3. Two Dimensional Example 



Consider the diagram of Figure 30, the following plots are the resulting mesh 



display of the decomposed coefficients three resolution levels down (m=-3) using a Haar 
wavelet and a constant phase-00 selection for each level. The top left diagram is the 
approximation coefficients, c,„. Going in a clockwise manner, the following plots equate 
to ^dn,, and ^d^, coefficients. 



Add one row of zeros for 
either: 




Phase-01 

Phase-11 




45 



wov/el«t Hor'iz: Maar wov'<5l«t Moor" 




O m 


D 


D2ro 


O 




Hor- iz p>Kos«: OOT \/er"t p> ^ 



R«solLJ'ti^m —3 



Figure 37. Mesh Display of the Coefficients for level m=-3 
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Figure 38. 



Contour Display of the Coefficients for level 
m=-3 
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Figures 39 and 40 compare the original and reconstructed image coefficients 



displayed as mesh and contour plots. 



Reconstructed C array Actual C array 




Resolution Level 0 



Figure 39. Reconstructed and Original Mesh Comparisons 



Reconstructed C array 




Actual C array 

100 1 ^ 



50 ■ 



50 100 



Resolution Level 0 

Figure 40. Reconstructed and Original Contour Comparisons 
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V. MULTIPLE PHASE DWT 



For both the one and two-dimensional transforms, we selected a single phase for 
decomposing to the next lower resolution level. If we decompose down to a fixed level 
m=M while picking different decomposition phases for each level, we generate a 
different set of coefficients. This set of decomposition coefficients and any other set of 
decomposition coefficients will successfully reconstruct the original input array as long 
as we note the phase used to obtain each particular resolution level from the previous 
one. Such a disparity in the sets of coefficients suggests that the DWT process of a two 
times decimation is not shift-invariant. In other words, if we decomposed the input array 
with one type of decomposition scheme, our coefficients would be different if the input 
array is shifted before the low pass and band pass filtering processes occurs. As a note, 
using the energy check routine, we can see the distribution of the energy of the 
coefficients change as the phase decomposition scheme varies. So this shift variant 
property can be further exploited by varying the phase for possibly optimizing or 
concentrating the energy of families of signals to certain resolution levels. 

However, we desire linear shift invariant (LSI) systems for the analysis of data. 
For the continuous WT, this property is generally true, however for the DWT case, this 
aspect is not. One way to bypass this obstacle is to account for all the phases during 
each decomposition. We call this technique the multiple-phase DWT (MP/DWT). A 
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thorough discussion for the one-dimensional case will first be presented, which will then 
be extended into the two-dimensional case. 

A. ONE-DIMENSIONAL MP/DWT 

Consider the case of going from the data input resolution m=0tom=-l. The 
coefficients that can be generated are one of two sets, depending on whether we selected 
the phase-0 or the phase- 1 case. We really need only one of the two sets to decompose 
to the next lower resolution level m=-2 to maintain perfect reconstruction. Also, we 
need to pick which phase to decompose from m=-l to m=-2. Thus for level m=-2, 
there are four possible phase combinations that can result: 00, 10, 01, and 11. The 
ordering of the subsequent phases is read from left to right to correspond to decomposing 
from level m=0 to m=-2 and we call this a phase vector. Thus [io]C. 2.5 would be 
interpreted as the set of approximation coefficients at resolution level minus two and 
decomposed with a phase vector of [10] with an index of five. In general, for any 
resolution level m < 0, the total number of phase combinations is 2 “. 

Looking again at going from the m=0 to m=-l case only, we can see that by 
sliding the filter coefficients by one instead of two units to the right, the resulting 
coefficients alternate on the phase selection starting with the first coefficient in the phase- 
0 set followed by the first coefficient in the phase- 1 set as shown on the next diagram. 
We determine the total number of coefficients in level m=-l just like a linear discrete 
convolution, which is the number of data points from level m=0 plus the length of the 
filter mask minus one ( | {cq} | -I- N - 1). 
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Figure 41. Alternating of Coefficients by Phase 

Going from level m=-l to m=-2 will require a slight modification before shifting 
the filter mask by one to the right as discussed in the previous paragraph. Since two 
sets of decomposition coefficients are interleaved, a zero must be appended between each 
coefficient in the filter mask for properly decomposing each set of coefficients. The first 
four values in the resulting set correspond to the phases; 00, 10, 01 and 11. This 
association repeats for each four values in level m=-2. In fact, the spacing of each 
coefficient in a desired phase is properly indexed by the same amount as the number of 
zeros padded for display purposes in the Chapter III.A.5. 

This observation can be generalized to any arbitrary level m. In this case, separate 
the filter coefficients by 2 “-l zeros before the conducting single shift to the right. By 
noting the resolution level m, the order of the phases can be determined as the binary bit 
reversal of the binary values counting from zero to 2 “-l. We show this characteristic 
as a phase-tree diagram for resolution levels m=0 to m=-3. 
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Resolution Level 




Figure 42. Phase-Tree Diagram for Resolution Levels m=0 to 
m=-3 



We now define new filter sequences, h„, and g„ (m <0) for decomposing the values 
c„u, to the next level Cn,., The revised filter sequences are defined as 



= [ H-N*2) (0(„ *(-W^1) (0(„ ... (0)„ HO) (0)„ Ml) 1 

(40) 

8 „, = [ ^(-^+2) {01 gi-N^l) ( 0 }„ ... iOl ^( 0 ) ( 0 }„ ^(1) ] 
where {0}„=2’"'-l length zeros vector 
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Thus, the MP/DWT equations are given as follows [Ref. 3: p. 3] [Ref. 4]; 



A 






( 41 ) 






and 



m-l,i 






( 42 ) 



t«-JV+2-(W-l)(2'"-l) 



The total number of coefficients also can be generalized in terms of the size of the 



input array, | do | , the size of the h or g mirror filters, N, and the current resolution 



level, m: 



|fj - kol * (W-l)(2-’-l) 



( 43 ) 



So enough memory must be allocated if we desire lower and lower resolutions, since the 
number of possible phases increases by a power of two each step downward. 

Consider a BPSK signal shown in Figure 43. Using a Haar wavelet and a phase 
vector of [11111111], a contour display of the energy in the approximation coefficients 
and of the detail coefficients results and are shown next. Notice that much of the 180° 
phase shifts are not readily detectable. 
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-n resolution be) 



BF=*SK Signa 




Figure 43. Sample BPSK Signal [Ref. 4] 



ZoorTned Contour Display of ths Appro xi mo tio n F^lnose: 1 1 1 1 1 




1 1 



fsj m to « r- of oontours: 1 O 

Figure 44. Approximation 
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r of DotoM 



P»»nos«: 11111111 




Figure 45. Detail |d^| 



Now for the MP/DWT case, the 180° phase shifts of the BPSK signal are now very 
apparent in both the approximation and detail coefficients: 



N/l cj 1 1 1 fD H o s o Contour Disploy of tHo 




Figure 46. MP/DWT 
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K/1 1 1 } |=> H o s e Cor-itOLjr Dl3|sla>/ of ttno D«toM 




Figure 47. MP/DWT |d^|^ 

Twice the energy of the coefficients for the current resolution level will equal the 
total energy on adjacent lower level’s approximation and detail coefficients, since there 
are two valid sets of coefficients for the reconstruction to the higher resolution. So the 
energy checking routine discussed in the one-dimensional DWT is still valid as long as 
we account for a factor of one half when going up one resolution level. 

The actual zero padding of the filter masks is not the most optimum method in the 
determination of the decomposition coefficients’ proper phase order. A much more 
efficient technique can be determined by noting the pattern that develops as Equations 
42 and 43 and writing it out explicitly. The following diagram shows the developing 
pattern for m=-2. 
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m 


0 0 




m = -2 


C(0) 


=a3c(0) 0 




L = |c 0 


c(l) 

c(2) 

c(3) 


=a3c(l) 0 

=a3c(2)4a2c(0) 

=a3c(3)4a2c(l) 


0 0 
0 0 
0 0 


)(L-1) 2 


'Z 

II 

II 


C(<3) 


=a3c(4)4a2c(2)+alc(0) 0 






C(5) 


=a3c(5)+a2c(3)+alc( 1 ) oj 






C(6) 


=a3c(6)+a2c(4)+alc(2)4a0c(0) 


where; 


C(7) 


=a3c(7)4a2c(5)+alc(3)+a0c(l) 


[aO al a2 a3 ... aN-1] 


C(L+K-3) 


=a3c(L- 1 ) +a2c(L- 


3)+alc(L-5 )4a0c(L-7) 


[h(-N+2) „.h(0) h(l)] 


C(L+N-2) 

C(L+N-1) 


= 0 a2c(L- 

= 0 a2c(L- 


2)+alc(L-4)+aDc(L-6) 
1 )+alc(L-3)+a0c(L-5) 


or 

[g(-N+2) g(0) g(l)] 


c(L^N) 


M 

O 

O 


alc(L-2)4a0c(L-4) 




c(L+2N-3) 


O 

O 

II 


alc(L-l )+a0c(L-3) 




C(L+2N-2) 


O 

O 

II 


0 


a0c(L-2) 




c(L+2N-l) 


= 0 0 


0 


aOc(L-l) 














(N-l) 









Figure 48. MP/DWT Pattern Development for m=-2 



Notice that the coefficient set, c„,.,, is a summation of N vectors with | Cq | +(N- 
1)(2 “-1) elements. We pad some of the vectors with zeros either in the beginning or the 
ending of the approximation sequence of the current level and multiply by an appropriate 
scalar value from one of the filter coefficients. The following equations summarizes the 
development; 



(«.-.) = E Op 

p^O 


Op = s/2hi-N+2+p) 


(44) 


N-l 


R-.} 

p=0 




(45) 
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where 



(46) 



K ‘ [{24 {«4 {a„)l 

{Onjb} = {zeros coefficient row vector of size (N-l-p)*2 ”} 

{0„e} = {zeros coefficient row vector of size p*2 '"}. 

This vectorized technique is faster than the iterative process of padding the filter 
coefficients appropriately. However, we need more buffer space to calculate the data, 
which is dependent on the number of resolution levels desired. 

The reconstruction process works in a very similar fashion as in the DWT case. 
All we need is the selection of the coefficient sets, which in turn determines the unique 
path back upwards in the phase-tree diagram, similar to Figure 8. Once we select the 
desired phase at the lowest resolution level, the we extract the corresponding 
approximation and detail coefficients and put them in a packed format. Then the DWT 
reconstruction commences to get up one level, where we extract the next set of detail 
coefficients. We repeat this process until we reach resolution level m=0. 

B. TWO-DIMENSIONAL MP/DWT 

Since we define the basis functions as separable in the two orthogonal directions, 
we can use the same arguments discussed for the one-dimensional case. When we talk 
about coefficients and masks, they are now two-dimensional arrays as discussed in 
Chapter IV. Thus, the basic equations for the two-dimensional MP/DWT are 



57 



y=-N**2-(N*-l)(2--l) 

*=-A/^*2-(A/,-l)(2-"-l) 



y.-iC".'’) = 2 5 ^ 5 ; Hp^(jjiyjfl*j-\.b*k-\) 

y=-W*+2-(A/*-l)(2--l) 

*=-A/,-2-(W^-l)(2'"-l) 



(<.,!>) =2 5 : 5 : ap,mcJfl*i-l,b*k-\) 

;=-N*+2-(A/»-l)(2--l) 

Jt--W^^2-(A^^-l)(2-"-l) 



’4., (a, 6). 2 5^5; G,G,(/^)c->y-l,6+t-l) 

/•=-A/*-2-(A/*-l)(2---l) 
t— A/,^2-(W^-l)(2-*-l) 



( 47 ) 



(48) 



(49) 



(50) 



where the top left comer of the mask corresponds to location (a,b)=(0,0), and the 
positive directions are to the right and downward [Ref. 3: pp. 2-3]:. 

Even the one-dimensional vectorized process converts into a "matricized" process. 
Here, the resulting coefficient array is the matrix addition of N^xN^ total matrices, each 
appropriately padded with zeros on all four sides and multiplied by one of NhXNy scalar 
values. The two-dimensional equations result as follows: 



N,-l Np\ 

^ ir ^pq 

p-0 q=0 



= 2H,H^(-N^2^p,-N^2^q) 



(51) 
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bp, = 2H,G^i-N^2^p,-N^2^q) 



(52) 



= E E 

p=0 ^=0 



N,-l N,-l 

{"4-,) = E E Z. = 2G,HJi-N*2*p,-N*2*g) «3) 

p-0 q-Q 

N,-l V> 

{’<?.->} = E E4. Z. /« ' 2GfiJ^-N*2*p,-N*2*q) 

p-0 q=0 

This time we define graphically in the following figure: 




Figure 49. Definition of Z„ 
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The phase-tree diagram can be viewed in three dimensions to more appropriately 
show the proper interleaving of the phases. Since there are four possible phases, there 
will be four more locations in the indexing of the data on the next lower resolution level. 
The total number of coefficients would be 

where RowSq and CoIsq are the number of row and columns of the original image. 

In order to view the phase-tree diagram correctly, we must consider each resolution level 
as a new plane for the ordering of the coefficients. 




m = 0 
m = -1 
m = -2 



Figure 50. Phase-Tree Diagram 
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The MP energy display of the coefficients closely resembles the display for the 
single phase case. This energy display differs by the fact that they are the normalized 
average energy values for each level. As before, we normalize the energy to the 
resolution level m=0. Consider as an example, a MP/DWT of the image in Chapter IV 
using Haar wavelets and decomposing down to m=-3. The following display shows the 
MP/DWT energy plot: 



wavelet Inariz: Hoar 




a 2 




mr-» uj 1 1 i p> 1 ^ a » e o a s e 



wavelet vert: Haar 




a 1 
a 3 



1 - 
o.s - 
o. e - 
o .^ - 
0.2 - 



Figure 51. MP/DWT Energy Display 



61 



and the energy check is as follows: 



Verification of the energy in each resolution level 



Level Energy in 

m level m-l= 

c+dl+d2+d3 



Energy in C Difference 

level m 



-2 0.441935911016949 

-1 0.751059322033898 

0 1 . 000000000000000 



0.441935911016949 

0.751059322033898 

1.000000000000000 



0.000000000000000 

0.000000000000000 

0.000000000000000 



The MP/DWT of the three-dimensional energy display for all the coefficients of 



the image and the corresponding contour images for the m=-3 case are as follows: 



wavelet horiz: Haar 




c 
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wavelet ven; Haar 



dl 

d3 






multiphase case 



Resolution level -3 



Figure 52. Mesh Display of the Stored Energy for m=-3 
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wavelet horiz: Haar 




d2 





d3 




Figure 53. Contour Display of the Stored Energy for m=-3 
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VI. MATLAB DWT ROUTINE DESCRIPTION 



The Appendix lists the developed DWT algorithms. We categorize them primarily 
into two classifications, the one and two-dimensional cases. In each situation, we further 
define two subgroups, the single phase DWT and the MP/DWT. Included with the 
software is an ASCII version of this chapter in a README.TXT file. This chapter first 
describes the system configuration to run the routines. Then we list brief steps for each 
of the two general cases. These steps are not totally inclusive, but the program’s 
prompts will fill any other gaps in the routine. Finally, we discuss important 
considerations for obtaining graphical output of the results since this is highly dependent 
on the system hardware. 

A. SYSTEM CONFIGURATION 

Minimum hardware requirements for the PC version are an Intel-386 
microprocessor with a math co-processor, Microsoft DOS or Digital Research DOS 5.0, 
VGA, at least three megabytes of hard disk space, and four megabytes of RAM, for most 
of the DWT applications. However, we recommend at least eight megabytes of RAM 
for using arrays greater than 10,000 elements, or if m<-5 in the two-dimensional 
MP/DWT case. Also, some DOS memory managers may conflict with the Pharlap DOS 
extenders that MATLAB Version 3.5 uses. For the Sun workstations, in order to operate 
PROMATLAB Version 3.5, we recommend the Open Windows or Sunview graphical 
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environments with at least eight megabytes of RAM. Any further memory can be 
negotiated with the system administrator. Note that UNIX-based operating systems 
discriminate characters between upper and lower case letters. 

B. GENERAL PROCEDURAL STEPS 

Initially, we must be in the MATLAB environment prior to running any of the 
routines. Either generate in MATLAB or load the sampled data for resolution level 
m=0, and assign a variable name. If and files with the prefix, "leg", and the suffix, 
".met" (all in lower-case letters), exists in the current directory, they will be deleted once 
the DWT routines commence. Thus if you desire to retain these files, rename them. 
The DWT routines use these files as the output graphic files for the routines. 

1. One-Dimensional General Procedures 

a. If you have no data, the routines include three files as data that can be loaded. 
Type "load" followed by a space and then one of the three names: "etbO.m", 
"et90.m", or "bpsk". The variable name will be either "et60", "et90", "bpsk". 

b. Invoke "wavld" in lower-case letters. 

c. Enter the variable name of the input data. 

d. Enter the sampling frequency, <RET> if not known. 

e. Select the desired h coefficients. 

f. Select "Y" if you desire the MP/DWT routine and skip to step k. 

g. In you did not select the MP/DWT routine, the single phase DWT commences, 
decomposing the coefficients until resolution level m=- riog 2 ( | c„ | )1 with an 
initial query of if you desire to see the intermediate resolution level displays. 

h. Choose either the phase-0 or phase- 1 decomposition for each resolution level. 
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i. After the final decomposition, the energy checks commence. 

j. Select "Y" if you desire to do the reconstruction process. If so, answer yes or no 
if you want to display the intermediate levels. 

k. Lastly, the time/ scale display commences, and the series of questions ask if you 
want to zoom in or out and what number of contour levels you desire. 

l. For the MP/DWT case, steps i to k are essentially the same except for the 
reconstruction, where you must initially select the desired phase, thereby fixing 
the phase path back up to level m=0. 

2. Two-Dimensional Procedures 

a. If you do not have any input data, invoke "idata" for a small selection of images. 
The variable "im" (in lower-case letters) will be the input for step c. 

b. Invoke "wav2d" in lower-case letters. 

c. Enter the variable name of the input data array. 

d. Select the desired h coefficients in the horizontal direction. 

e. Select the desired h coefficients in the vertical direction. 

f. Select "Y" if you desire the MP/DWT routine. 

g. Enter the number of resolution levels to go down. An initial good selection is 
four since this will cover all four possible phases. 

h. For the single phase DWT case, select one of four of the decomposition phases 
for each level. 

i. After the final decomposition, the energy checks commence. 

j. For the single phase DWT case, select "Y" if you desire the reconstruction 
routine. 
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C. OBTAINING GRAPHICAL OUTPUT 



Except for the initial time/frequency plots, all other displays will query you if you 
want to store the plot. If so, the plot is stored as a metafile with the prefix of "leg" and 
the suffix ".met". A number starting from zero is put between the prefix and suffix, for 
example, leg 14. met would correspond to the fifteenth stored plot in the DWT routine. 
Output of these metafiles are hardware and system dependent. Refer to the GPP 
command in the MATLAB User’s Guide for more detailed information. Two 
C-shell script files, called "metal3" or "metal4", will generate a temporary Postscript 
file for plotting on a Postscript printer all the metafiles retained from the DWT routine. 
The use of these files is currently only guaranteed to work at the US Naval Postgraduate 
School Electrical Engineering Department’s Sun workstations, and we list them here as 
a guide in generating other C-shell script files for a particular system. For the DOS 
version, the following command will generate all the metafiles onto a HP LaserJet 111 
printer as long as the we invoke the following DOS command in the same directory as 
the metafiles: for %f in (leg*. met) do gpp386 %f /djetl50 /ol /fpm . If you put it in 
a DOS batch file, echo the percent sign ( %% instead of % ). Please note that the 
computer may take some time generating these plots. Be also aware that depending on 
the size of the plot, the graphic files may exceed the printer buffer memory. 
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VII. CONCLUSIONS 



In this thesis, we developed discrete wavelet transform algorithms for both the one 
and two-dimensional cases. The iterative mechanism of decomposing the data 
coefficients by a convolution-and-subsample process is achievable. These realizations 
require that the scaling and wavelet functions form orthonormal sets with compact 
support. Also for the two-dimensional case, the basis functions were separable in the 
two orthogonal directions. 

We found that the output of the decomposition coefficients can be causal provided 
that the definitions of the indices for both filter coefficients, h and g, must be selected 
to accommodate the causal condition. In fact, the definition of the g filter coefficients 
must be adjusted from the definitions given by Mallat and Daubechies. We found for 
the one-dimensional case, the decomposition is simply a dot product of the filter 
coefficients with a set of values windowed from the higher resolution level approximation 
coefficients. This window subsequently shifts two units to the right for the next 
decomposed coefficient. We found that the reconstruction process uses the same 
decomposition algorithm if we separate each of the individual input values by a zero, the 
filter coefficients are time- reversed, and the shift is now one unit to the right at a time. 

The two-dimensional case was a natural extension of the one-dimensional process. 
We found that the decomposition can be realized as a mask element multiplication and 
summation, with the mask shifting two units to the right until the row completion and 
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then two units down to process the next row of coefficients. The shifting of the mask 
is similar to a raster scan. The reconstruction follows the same process as in the 
decomposition but each of the input array coefficients are separated by a row and a 
column of zeros, the filter masks are spatially reversed in both direction, and the raster 
shift is one unit instead of two. 

We found that there were two possible phases in decomposing the data in the one- 
dimensional case, and four possible phases for the two-dimensional case. This led to the 
conclusion the current decomposition algorithm was very shift variant, contrary to the 
continuous version to the wavelet transform, which is shift invariant. We then 
determined that the discrete case can possibly approximate the continuous version if we 
account for all the phases during the decomposition. We then developed the multiple- 
phase discrete wavelet transform for both one and two dimensions. The coefficient sets 
for a particular phase can be interleaved with other sets of a different phase. We found 
that the multiple-phase was more capable in detecting discontinuous properties of data 
than in the regular discrete wavelet transform, such as a binary phase shifting function. 



69 



APPENDIX 



DWT MATLAB HLE LISTINGS 



A. ONE-DIMENSIONAL ROUTINES 



% IS SEP 92 

% WAVID.M The One Dimensional Discrete Wavelet Transform 

% 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: leg asp i@ece.nps. navy, mil 

% lam@cce.np8.navy.mil 

% Phone: (408) 646-3044/2772 

% 

% 

% This routine docs a one-dimensional dicrcte wavelet decomposition 
% of a one-dimensional data array. The algorithm allows for the 
% selection of the desired phase for the respective dccompositon 
% level. The following routines arc necessary for the proper 
% operation of this algorithm: 

% 

% chkputcr.m — algorithm checks for compatible hardware 
% daubdata.m — determines the h coefficients for a daubechics 
% compactly supported orthogonal wavelet 

% mp.m — Multiphase decomposition 
% phsl.m — selects the phase-1 decomposition 

% phsO.m — selects the phase-0 decomposition 

% dpscocf.m — display the coefficients for each level 
% cnrgld.m — determines the energy check for each resolution 
% rcconld.m — reconstruction algorithm for verification 
% plane. m — displays the coefncienls in the time-scale domain 
% strplt.m — stores the plots in meta-Hles 

% 

clc 

chkputer % Checking for the basic hardware necessary 



% 

% cO corresponds to the resolution level 0, the original data 

% array. 

clc 

c0 = input(’ Enter the name of the data in vector row format: ’); 

% We will check if the data is in row format, if in column format, 
% take the transpose of the input. 



(il jl]=size(c0); 

if jl = = l, c0 = c0*; 
elseif il > 1 & jl > 1, 

disp(*Bad Input Data, Restart the Program’) 
return 



end 

clear il jl 
% 
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f8amp=input(’ Enter the sampling frequency (Hz): ’); 



% 

% The following code determines the desired wavelet 

for i = 1 :8,dUp([’ *]),cnd 

ql =[’Enlcr the number for the desired choice: ’ 

’(1) Haar h coefficients ’ 

’(2) Daubechies h coefficients ’ 

’(3) Own set of h coefficients ’ 

1 ; 

disp(ql) 
qqq = input(’ ’); 

if qqq==3. 

h = input(’ Enter your own set of h coefficients’); 

%Chccking if the data are compact support 
a = sum(h);b=h*b’; 
to!l = lc-12; 

while abs(l-a) >toll & abs(.5-b) > toll 
disp(’Your h coefficients are NOT compactly supported!’) 
h = input(’ Re-enter the h coefficients or ctrl-z to stop’); 
end 

wwsvlet = ’OTHER’; 
elseif qqq==2, 

qq = input(’ How many taps (2-10)? ’); 
h=daubdata(qq); 

wwavelel = [’Daubechies ’ ,num28ti(qq) , ’-tap ’] ; 
else h = [.5 .5];wwavelet = ’HAAR’; 
end 



h = 8qit(2)*h; % The factor of sqn(2) normalizes the energy 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 



% 

% We will intialize some of the constants: 

LL=length(cO); 

Nh = length(h); 

lowest = -ccil(log(LL)/log(2)); 

lvl = 0; 

pllcnt=0; 

% 



Some definitions of variables used throughout the entire 
algorithm: 

wwavelet — name of the selected wavelet 
LL - the length of the input array 
Nh — the number of h and thus g coefficients 
lowest — the number of resolution levels below the input 
pltcnt — plot counter for storing desired plots 
h — the "h" coefficients 

g — the *g* coefficients 

lowest — the lowest resolution level 

phsvct - the phase vector for recording the appropriate phase 
zro — zero vector to record the appropriate zero padding 
for each resolution level. 
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% 

% Generating the g coefficients 

g = fliplr(h); 
for i=2:2;Nh 
g(l,i) = -gCl.i): 
end 

% 



% 

% The following set of Lines branches to the multiphase 
% decomposition algorithm if desired. Program ends after 
% the multiphase execution 
clc 

q = input('Doyou desire the muliphase decomposition (Y/N)? *.*8*); 
ifq==’Y' I q==y 
mp 
return 
end 

% 



% 

% The Decomposition Routine 

% Determine if we want lo ignore the display for each resolution 
% level’s coefficients. 

pitl =input( ’Bypass the display of the Coefficients (Y/N)? ’.’s’); 
clc 

phs=2; % This just sets "phs" initially lo be NOT equal to 1 or 0 
phsvet = [] ; 
zro = 0; 



% 

% The following code sets the number of variables necessary 
% to record the coefficients, ie. lowcst=-2, we have 
% c0,cl.c2,dl,d2 

LLL=LL; 

for lvl = -l:-l:lowcst 
while 1 = = 1 

phs = input(’ Enter the desired Phase [0/1] for this resolution level: ’); 
ifphs= = l I phs = =0, break, end 
end 

phsvci = [phs phsvcl] ; 

eval( [ ’ c ’ ,num2str(abs(lv 1)) , ’ = phs ’ ,num2str(phs) , ’( c ' , . . . 
num2slr(abs(lvl)- 1) , ’,h) ; ’]) 

e val( [ ’d ’ ,num2slr(abs( Ivl)) , ’ = phs * , num2str(phs) , ’(c ’ , . . . 
num2str(abs(lvl)- 1), ’ ,g) ; ’]) 

% This if statement runs the display of the coefficients if the flag 

% pltl is set to NO 

ifplll = = ’N’ I plll = = ’n’, dspeoef, end 

phs=2; 

eval([’LLL=max(LLL,length(c’ ,num2str(-lvl) , ’) •2‘’(-lvl)); *]) 
zro = [zro 2*(-lvl)-l]; 
end 

% 

% We now do an energy check of the coefficients 
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enrgid 

clg 

% 

% 

% We now run Ihc reconstruction option 

recon Id 

% 



% 

% We now display the coefficients in the time-scale domain 

plane 

% 

% END OF WAVlD.M 



% 

% Phase Zero decomposition for the 1-D case 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail; legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

% 

% phsO.m is a function m-file the does a phase 0 decomposition of 
% of the data array using either the h or g coefficients. The main 
% program that calls this routine is "wavld.m" 

% 

% The input arguments for this routine: 

% 

% data — input data 

% h — the h or g coefficients 

% L — the length of the data vector 

% H — the length of the h or g coefficient vector 

% 



function x=ph80(data,b) 

L=length(data); 

H = length(h); 



% The following 4 lines check if the data array is even, for the 
% phase 0 case, we must pad it with one zero if true 
if rem(L/2,floor(L/2)) = = 0 
data = [data 0]; 

L=L+1; 

end 

d ata = [zeros( l,H-l)data zeros( 1 ,H-2) ] ; 

L=L+2*H-3; 
a = 0; 

for ii = l;2:Lr(H-l) 
a = a+l; 

x(l,a)=data(ii:ii + H-l)*h’; 
end 
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% 

% Phase One decomposition for the 1>D case 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail; legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408)646-3044/2772 

% 

% phsl .m is a function m-file the does a phase 1 decomposition of 
% of the data array using either the h or g coefficients. The main 
% program that calls this routine is "wavld.m" 

% 

% The input arguments for this routine: 

% 

% data - input data 

% h " the h or g coefficients 

% L — the length of the data vector 

% H - the length of the h or g coefficient vector 

% 

function x=phsl(data,h) 

L=length(data); 

H = length(h); 

% The following four lines checks if the data vector is odd. if it 
% true for the phase 1 case, we must pad it with one zero 
if rem(172,floor(I72)) ~ = 0 
data = [data 0]; 

L=L+1; 

end 

data = [zeros(l,H-2) data zeros(l ,H-2)]; 

L=L+2*H-4; 

a = 0; 

for ii=l:2:L-(H-l) 
a=a-M; 

x(l,a)=data(l,ii:ii+H-l)*h’; 

end 



% 

% dspcoef.m Display the coefficients 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School. Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone; (408)646-3044/2772 

% 

% dspcoef.m dispUys the coefficients for each particular 
% resolution level with the appropriate scaling factor. There is 
% also a special routine embedded in this algorithm for the 
% proper bar display of the HAAR case. The main program is "wavld.m’ 
% 

% The variables arc defined as follows: 

% 

% z - the number of zeros to be padded blwn coeff s 

% Ivl — variable from *w.m* for resolution level 
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% fsamp -* sampling frequency of the data 
% wwavelet — string for the selected wavelet 



% 

clg 

hold off 
axis(’ normal’) 

% For the HAAR case do the following: 
if wwavclcl(l) == ’H’ j wwavclet(l) == ’h’ 
n = 2'^(-lvl); 

cval( [ ’ctcmp = 2 Iv 1/2) *c ’ , num2sir(-lv 1) , ’ ; ’)) 
cval([’dtemp = 2"" (lvl/2) *d ’ ,num2str(-lvl) ,’;’]) 
ymax = max(max(clcmp) , max(dtemp)) ; 
if ymax <0, ymax = 0;end 
ymin = min(min(ctemp) , min(dlemp)); 
if ymin>0,ymin = 0;cnd 
L= length(ctcmp) ; 
if L= = l 

axis([0 (n* 1.5-1) /fsamp 1.2*ymin 1.2*ymaxJ); 
subplot(21 1) 

plot([0 0 n* 1.5-1 n*1.5-l]/fsamp,[0ctemp ctemp 0]) 
lille([ ’Approximation at Resolution level ’,num2slr(lvl)J) 
xlabcK’Haar wavelet’) 

axis([0 (n*1.5-l)/fsamp 1.2*ymin 1.2*ymax]); 
subplot(212) 

plot([0 0 n* 1.5-1 n*1.5-l)/fsamp,[0dtemp dtemp 0]) 
Utle([’Dctail Phase- ’,num2slr(phs)]) 

pause 
strplt 

else 

x = [0;n:n*L-lJ + .5*n; 
x=x/fsamp; 

axis([0 x(length(x)) 1.2*ymin 1.2*ymax]); 
subplot(21 l),bar(x, ctemp) 

UUe([ ’Approximation at Resolution level ’,num2slr(lvl))) 
xlabeK’Haar Wavelet’) 
axis([0 x(length(x)) 1.2*ymin 1.2*ymax)); 
subplol(2 1 2) ,bar(x .dtemp) 

tiUe([’ Detail Phase- ’,num2sir(phs))) 

pause 
strplt 
end 

% For NON-HAAR cases, do the following: 
else 

z = 2>lvl)-l; 

eval([ ’ctemp = 2*(lvl/2) *ptzero(c’ ,num2slr(-lvl), ’ ,z) ; ’]) 

cval([’dtcmp =2'‘(lvl/2)*ptzero(d’ ,num2str(-lvl), ’ ,z); ’]) 

ymax = max(max(ctemp),max(d temp)); 

if ymax < 0, ymax = 0; end 

ymin = min(min(ctemp) , min(dtemp)) ; 

if ymin> 0, ymin =0; end 

% We can display the data with a window size of the original 
% input, or show the spillover values 
q = [’Enter the type of display:’ 

’ 1 - Windowed Display ’ 

’ 2 - Non-Windowed ’]; 

clc 

disp(q) 
dUp([’ ’]’) 
q = input(’ ’); 
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if q==2 

X = [0: length(clcmp)- 1 ] + . 5 ; 
x = x/f8amp; 

axls([0 x(lenglh(x)) ymin ymax]); 
subplot(2 11), bar(x .cterap) 

liUe([ ’Approximation al Resolution level ’,num28tr(lvl)]) 
xlabel([wwavelet,’ Wavelet’]) 
axi8([0 x(length(x)) ymin ymax]); 

8ubplot(2 1 2) ,bar(x,dtcmp) 
tiUe([’ Detail Phase- ’,num28lr(phs)]) 

xlabel([’Non-Windowcd display -- max windowed value; 
num28tr(length(x)- 1)]) 

pause 

strplt 

else 

x = [0:LH] + .5; 
x^x/fsamp; 

axis([0 x(length(x)) ymin ymax]); 
subplot(21 l),bar(x,ctcmp(l,l:LL)) 
tiUe([ ’Approximation at Resolution level *,num28tr(lvl)]) 
xlabel([wwavclet,’ Wavelet’]) 
axis([0 x(length(x)) ymin ymax]); 
subplot(2 1 2).bar(x,dlemp( 1,1: IX)) 
titled ’Detail Phase- ’,num28lr<phs)]) 

xlabel([’ Windowed display -- max windowed value; ’.num2slr(lX-l)]) 
pause 
stxplt 
end 



^••••*******************«***********«******************************* 

% enrgld.m Detennining the Energy in each resolution level 
% One Dimensional Case 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: lcgaspi@cce.np8.navy.mil 

% lam@ecc.np8.navy.mil 

% Phone: (408) 646-3044/2772 

% 

% For the aingle phase case, the energy is normalized to resolution 
% level 0. All the energies of the ’’c’l” and the "d’s” of a 
% particular level should add up to the energy of the "c’s" of the 
% next higher level. 

% 

% enrge - energy array of the "c’s" 

% enrgd — energy array of the "d’s" 

% norm — energy of the data "cO" 

% esura — energy sum array (enrge + enrgd) 

% 

% The calling routine is ’wavld.m" . 

% strplt. m is a necessary m-filc 
% 

enrgc =zcro8( 1 , -lowest) ; 
enrgd =zeros(l , -lowest); 
a=0; 

norm=8um(c0 .^2); 
for ii = lowest: 1:-1 
a = a+ 1; 

eval([ ’cnrgc(a) =8um(c’ ,num2str(-ii) , '^2) ; ’]) 
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e va 1( [ ’ cn rgd (a) = 8Uin(d ’ , num28lr(- ii) , ’ . 2) ; ’ ] ) 
end 

enrgc = cnrgc/norm; 
enrgc = [enrgc 1]; 
enrgd = cnrgd/norm; 
cig 

axis((lowc8l-2 1 0 1.2]); 

8ubpIot(2 1 1 ) ,bar([Iowe«t: 0] . [enrgc] ) 
liUe([ ’Normalized energy of the *c’ coefficients Phase: 
8elstr(fliplr(phs vcl + 4 8))] ) 
xIabeIC Resolution level’) 
axis([loweat-2 1 0 1.2]); 

8ubplol(2 1 2), bar<[ lowest:-!], [enrgd]) 

title(’Nonnali 2 ed energy of the ’d’ coefficients’) 

xlabel(’ Resolution level’) 

pause 

strplt 

axis 

% Energy check of the previous plot 

esum = enrgc + [enrgd 0] ; 
cic 

dd = [’ Verification of the energy in each resolution level 

’Level Energy in Energy in C Difference’ 

’ level m-1 level m 

’ c+d 

• ]; 

disp(dd) 
j = -lowest; 
for i= 1: -lowest 

a = esum(i) ; b = en rgc( 1 , i + 1 ) ;c = 1 -j ;d = abs(a-b) ; 
fprintf(’%2.0f %17.15f %l7.15f ’.c.a.b) 

fprimf( ’ % 1 7. 1 5 fin ’ ,d) 

end 

pause 

cIc 

axis; 



% 

% reconld.m One Dimensional Single Phase Reconsmetion 

% 



% By: 
% 

% 

% 

% 

% 

% 



LT J. E. Legaspi 
Professor Lam 

US Naval Postgraduate School, Monterey CA 
e-mail: legaspi@ece.nps.navy.mil 
lam@ece.nps.navy.mil 
Phone: (408) 646-3044/2772 



% This routine conducts a reconstruction of the the data with the 
% proper phase selected between resolution levels. This is a sub- 
% routine for the main wavelet program ’wavld.m’. 

% 

% The following variables are further defined: 

% 

-- reconstruction coefficients 

— defined in ’wavld.m’ 

— reconstructed c’values at a particular Ivl 

— phase vector defined in "wavld.m". 



hh & gg 
lowest 
cwork 
phsvcl 



% 

% The following m-filcs arc ncccsaary for proper operation: 
% 

% wavld.m — main program 

% rphfll.m — rcconslruction phase- 1 function m-filc 

% rphsO.m — reconstruction phasc-0 function m-file 

% strpll.m — storage of the displays 

% 



clc; 

disp([’ ’]’) 

q = input(’Do you desire to do the Reconstruction Comparison (y/N)?’,’s’); 
clc 

ifq== *Y‘ I q== y. 

% generation of the reconstruction coefficients 
hh = fliplr(h) ; gg = fUplr<g) ; 
u = l; 

eval([’c work = c‘ ,num2sli<- lowest), ’ ; ’]) 
for lvl = lowesl:l:-l 
eval([’dwork=d’,num2str(-lvl),’;’]) 

cval([’cwork=rphs*,num2str(phsvct(ii)),’(cwork,dwork,hh,gg);’J) 
eval([ *a = lcngth(c’ ,num2sli(-lvl- 1 ), ’); ’]) 
cwo rk = cwork( 1 : a) ; 
ii=ii+l; 

phswrk = phsv ct( 1 : leng th(phsvct)- 1 ) ; 

% Comparison of the Original and the Reconstructed data 
clg 

subplot(21 1) 

eval([’bar<c’,num2slr<-lvl-l),’)’J) 
titJc([’ Original Data Phase: ’,sctsti<(phswrk+48))]) 
xlabel([* Wavelet; \wwavelet]) 
subploK212) 
bar(cwork); 

tiUcC Reconstructed Data’) 
xlabel([’Rcsolution level ’,num2str<lvH- 1)]) 
pause 
strplt 
clg 

a = eval([’c’,num2sU<-lvl-l),’-cwork’]); 
bar(abs(a)); 

tiUe( ’Absolute Error in the Reconstruction’) 

xlabel([’Resolution Level ’,num2str(lvl+ 1)]) 

pause 

strplt 

clg 

phswrk = phsw rk( 1 : lengtfa(phs wrk)- 1 ) ; 
end 
end 



function oul = rphsO(cwork,dwork,hh,gg) 

% rphsO.m Reconstruction Phase Zero routine 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.np8.navy.mil 

% Phone: (408) 646-3044/2772 

^***********************«****«*«*«*«, •«•*••« 
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% 

% This routine does a reconstruction of a particular level with a 
% phase 0 selection. The calling program is *wavld.m'. Input 
% arguments are: 

% 

% cwork — ’c’ coefficients of a particular level 

% dwork — "d* coefficients of a particular level 

% hh — reconstruction coeff 8 for the ’c’s" 

% gg " reconatrucilon coeffs for the *d’s* 

% 

% "out" is the output argumemt that equals the c coeffs of the 
% next higher resolution level. 

% 

L= length(cwork) ; 

xl=D; 

x2=0; 

for j = l:L 

xl =[xl cwork(l,j) 0]; 
x2 = [x2dwork(l,j)0]; 
end 

xl=[xl zcros(l ,length(hh)-2)]; 
x2 = [x2 zeros(l,length(gg)-2)]; 
a = 0; 

for j = 1 :length(xl)-(length(hh)-l) 
a = a+ 1 ; 

out( 1 ,a) = xl ( 1 ,j:j + length(hh)- 1 ) *hh ’ + x2( 1 j:j + length(gg)- 1 ) *g g ’ ; 
end 



function out = rphsl(cwork,dwork,hh,gg) 

^•*****************************,******, *,**•••,,*•*, 

% rphsl.m Reconstruction Phase One routine 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi(^ece. nps.navy.mil 

% lam@cce.nps.navy.mil 

% Phone: (408) 646-3044/2772 

% * 

% 

% This routine does a reconstruction of a particular level with a 
% phase 1 selection. The calling program is "wavld.m". Input 
% arguments are: 

% 

% cwork — *c" coefficients of a particular level 

% dwork " *d" coefficients of a particular level 

% hh " reconstruction coeff s for the "c’s* 

% gg " reconstruciton coeff s for the "d’s" 

% 

% "out* is the output argumemt that equals the c coeff s of the 
% next higher resolution level. 

% 

L=length(cwork); 
xl=D; 
x2 = U; 
for j = l:L 

xl =[xl cwork(l,j) 0]; 
x2=[x2 dwork(l ,j) 0]; 
end 

xl =[0 xl zeros(l .length (hh) -2)]; 
x2 = [0 x2 zcros(l ,lcngth(gg)-2)]; 
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a=0; 

for j = 1 . lcngth(xl)-(lcngth(hh)- 1) 
a=a+ 1 ; 

oul( 1 ,a) = xl( 1 J:j + length(hh)- l)*hh' + x2( 1 ,j:j + length(gg)- l)*gg ’ ; 
end 



* 

% plane. m develop the phase plane diagram 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School. Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam@ecc. nps.navy.mil 

% Phone: (408) 646-3044/2772 

^••*****»***«*****»*»*********************************************** 

% Phase Plane Diagram for the One Dimensional Single Hiase case 
% This routine is called by the parent routine "wavld.m”. All the 
% variables defmed in "wavld.m* are used here. New variables: 

% 

% c " (1 -lowest) X LL ’c’ coefficient matrix 

% d — (-lowest) X LL "d" coefficient matrix 
% N,N1 -- variables for the number of contour levels 
% 

% The main subroutine used is ptzero.m funtion m-file 

% 

clc 

clg 

% Set up the c and d matrices from the "c*" and "d*" variables 

c = zeros(- lowest 1 , LL) ; 
d =zeros(-lowest,LL); 
c(l.l:length(c0)) = c0; 
clear cwork dwork 
for lvl = l:-lowc8t 

cwork = eval(( ’ptzero(c’ ,num2str(lvl) . ' .zro(lvl)) ']); 
dwork = cval([ 'ptzero(d ’ .num2str(lvl), ’ ,zro(lvl)) ’]); 
c(lv 1 1 , 1 : length(cwork)) = cwork ; 
d(lvl, 1 : lenglh(dwork)) =d work; 
end 

c = flipud(c.''2); 
d = flipud(d.^2); 

(i j]=size(c); 

x=0:j-l;yc = 0:i-l;yd = l:i-l; 
mesh(c); 

titled 'Energy Display of Approximation Phase: ',... 
selstr(fliplr(phsvct -♦- 4 8))] ) ; 

pause 
N = 10; 

contouffc.lO.x.yc); 

titled ’Contour Display of the Approximation Phase; ’,... 

sets tr(fliplr(phs vet -1- 4 8))] ) ; 
ylabel(’-n resolution level’); 
xlabel(’10 contour levels’) 
pause 
clg 

ql = inpul(’ Do you desire to zoom in on the display (Y/N)?’,’s’); 
while ql = = ’Y’ | ql = = ’y’, 
q2 = input(’Doyou want to see the plots again (Y/N)?’,’s’); 
ifq2==’Y’ I q2 == ’y’. 

mesh(c);title(’Previou8 Energy Display of the Approximation’); 
pause 
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conlour(c,N,x,yc);tille(’PreviouflContour Display of the Approximation’); 
ylabel(’-n resolution level’); 
xlabel([’Numberof coumours: ’,num28lr(N)l) 
pause 
end 
clc 

disp([’ *]’) 

disp([’x range: 0 to ’,num2slr(j-l)]) 
disp([’y range: 0 to -’,num28tr(i-l)]);di8p([ ]’) 

xl =input(’ Enter the minimum x-value;*) + 1 ; 
x2 = input(’ Enter the maxumum x-value:’) + l; 

yl =input(’ Enter the minimum *magnitude* y-value (least negative): ’)+ 1 ; 

y2 = input(’ Enter the maxumum "magnitude" y-value (mo8t negative) :’) + 1 ; 

y 1 = abs(y 1 ) ;y 2 = abs(y 2) ; 

y = l:l:i;y = fliplr(y); 

yll=y(y2); 

y22=y(yl); 

clg 

mcsh(c(y 1 1 : y 22 , x 1 : x2)) 

tiUe([’ Zoomed Energy Display of the Approximation Phase: 

8el8tr( fliplr(phsvct + 48))]) ; 

pause 

strplt 

clg 

Nl = 10; 
while N 1 < 999 
N = N1; 

contour(c(y n :y22,xl ;x2),Nl ,xl :x2,y 1-1 :y2*l) 
tiUe([ ’Zoomed Contour Display of the Approximation Phase: 
sel8tr( fliplr(phsvct + 4 8))]) ; 
ylabcl(’-n resolution level’) 
xlabel([’Numberof contours: ’,num28tr(Nl)]) 
pause 
suplt 

N1 =input(’ Enter the number of contour leveb (999 to continue) ’); 
end 

ql =input(’Zoomin or out Further (Y/N)? ’.’s’); 
end 
clc 

% Now the Detail Signal 
mesh(d); 

tille([’ Energy Display of the Detail Phase: ’,set8tr(fliplr(ph8vct+48))]); 

pause 

clg 

contour(d,10,x,yd); 

titie([*Contour Display of the Detail Phase: ’,setstr(fliplr(phsvct+48))]); 

ylabcl(’-n resolution level’); 

xlabel(*lO contour levels’) 

pause 

clg 

clc 

ql =input(*Doyou desire to zoom in on the display (Y/N)?’, ’s’); 
while ql = =’Y’ 1 ql = = ’y’. 

q2 = input(*Doyou want to see the plots again (Y/N)?’, ’s’); 
ifq2==’Y’ t q2 ==’y’. 

mcah(d);tille(’ Previous Energy Display of the Detail’); 
pause 

contour(d,N,x,yd);titic(’PreviousConlour Display of the Detail’); 
ylabel(’-n resolution level'); 
xlabel([’Number of contours: ’,num2str(N)J) 
pause 
end 
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disp([’x range: 0 to \num28lr(j-l)]) 
disp([’y range: *1 to *\num28U^i-l)]);di8p([ ]’) 

xl =inpul(’ Enter the minimum x-value:’)+ 1 i 
x2 = input(’ Enter the maxumum x*value:’) + l; 

yl =input(’ Enter the minimum ’magnitude" y-value (least negative):’); 
if yl = = 0, yl = 1; end; 

y2 = inpul(’Enterthe maxumum "magnitude* y-value (most negative):’); 

yl =abs(yl);y2 = ab8(y2); 

y = l:l:i-l;y = fliplr(y); 

yll=y(y2); 

y22=y(yl); 

clg 

meah(d(yl I:y22,xl:x2)) 

title([’ Zoomed Energy Display of the Detail Phase: 
8el8lf(fliplr(phsvct+ 48))]); 

pause 

suplt 

clg 

while N <999 

contour(d(yl l:y22,xl:x2),N,xl:x2,yl:y2) 
title([’ Zoomed Contour Display of the Detail Phase: 
8el8lj<fliplr(ph8vct + 4 8))]); 
ylabel(’-n resolution level’) 
xlabel([ ’Number of contours: ’,num28tr(N)]) 
pause 
suplt 

N = input( ’Enter the number of contour leveb (999 to end): ’); 
end 

ql =input(’Zoomin or out further (Y/N)? ’.’s’); 



% 

% Multiphase routine for The Discrete Wavelet Decomposition Routine 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

% 

% 

% This routine does a multiphase discrete wavelet decomposition of 
% a one-dimensional data vector. This routine requires the following 
% M-files; 

% wavld.m 

% sUplt.m 

% zOpad.m 



maxSamp =(Nh- 1 )*(2 '"(-lowest)- 1 )+ LL; 
coeffc=zero8(-lowest-H 1 .maxSamp); 
coeffc(-lowest-H 1 ,1 :LL)=cO; 
coeffd =zeros(-lo west, maxSamp); 
numcocf =zeros(- lowest + 1,1); 

% 

% 

% The Main Multi-Phase Decomposition Algorithm *** 
numcoef(- lowest -H 1) = LL; 

i = l; 

for row = -lowest:-l :1 

numcocf(row) = (Nh-1 )*(2“ i- 1 ) -f- LL; 
pi =zOpad(coeffc(row-H,l:numcoef(row-Hl)).h,i); 
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cocffc(row, 1 : length(pl)) = pl ; 

p2=zOpad(cocffc(row + 1 ,l:numcocf(row+ l)),g,i); 
cocffd(row,l:lcngth(p2)) =p2; 
clear pi p2 
i = i+l; 
end 

% 

% Display of Coefficients at each resolution level 
clc 

q = input(’Do you desire to sec values at d iff resolution levels (Y/N)? ’,’s’); 
if q = = ’Y’ I q = = ’y' 
plot([0:LH],c0/*') 

xlabcl([’n sampling points n*’,num2str(l/f8amp),’ seconds’]); 

ylabcl( ’Amplitude*) 

tiUe(’ Original Data (resolution level 0) Multiphase case’); 

pause 

strplt 

clg 

index=[-lowcst;*l:l]; 
xl=0;x2 = -lowcst+l; 
clc 

disp([ ’Levels to pick are from -1 to ’,num2str( lowest)]) 
while xl < 1 I x2 > -lowest 

xl =input( ’Enter first resolution level of the desired range: ’); 
x2 = input(’Enter second resolution level of the desired range: ’); 
xl=abs(xl);x2=abs(x2); 
if xl >x2, tcmp = x2; x2 = xl; xl =lemp; end 
end 

while xl< =x2 

pts = numcoef(index(xl)); 

Tc=coeffc(index(xl), 1 :pts); 

Td = coeffd( index( X 1 ) , 1 : pts) ; 

minpts =min(nun(Tc),inin(Td));ifminpt8 > O.minpts =0;end 
maxpts = max(max(Tc), ma x(Td)) ; ifmaxpts < 0 .maxpts = 0 ; end 
if minpts = = 0 & maxpts = =0, maxpts = .5 ;cnd 
u = [0pts 1.2*minpts 1.2*maxpts]; 
axis(u); 

subplot(2 1 l).plot([0:pts- 1 ] ,cocfTc(indcx(xl ) , 1 ; pts) , ’ • ’) 

titlc([’Multiphase Approximation Coeff at Resolution Level -’,num2str(xl)]) 

axis(u); 

subplot(2 1 2) ,plot([0 : pts- 1 ] ,cocffd(index( X 1 ) , 1 : pts) , ’ • ’ ) 
title(’Multipha8e Detail Coefficients’) 

xlabel([’n sampling points n*’,num2str(l/fsamp), ’seconds’]); 

pause 

strplt 

clg;clc 

axis; 

xl=xl + l; 
end 
end 

% 

% Multi-Phase Energy Check 

enrgldmp 

% 

% Phase Plane Determination 

c = cocffc.'^2; 

d=coeffd.^2; 

[i j] =sizc(cocffc); 
x=0:U^l; 
yc = 0:i-l ; 
yd = l:i-l; 

cc = c( 1 : - lowest + 1,1; LL) ; 
clc 

mcsh(cc); 

titlcC Multiphase Energy Display of the Approximation’) 

pause 

strplt 
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clg 

N = 10; 

conlouKcc. 10,x,yc); 

lille('Mulliphase Contour Display of the Approximation’); 

ylabel(’-n resolution level’); 

xlabeK’lO contour levels’) 

pause 

strplt 

clg 

ql =input(’Doyou desire to zoom in/out on the display (y/N)?’.’s’); 
while ql = =’Y’ | ql == ’y’, 

q2 = input(’Doyou want to see the plots again (y/N)?’,’s’); 
ifq2==’r t q2 == ’y’. 
mesh(cc) 

liile(’ Previous Multiphase Energy Display of the Approximation’); 
pause 

conlour(cc , N , x ,y c) 

tiile(’ Previous Multiphase Contour Display of the Approximation*); 
ylabel(’-n resolution level*); 
xlabel([ ’Number of contours; ’,num2str(N))) 
pause 
end 
clc 

di8p([* *]*) 

dispd’x r^ge: 0 to *,num2str(j-l)]) 

dispd’y range: 0 to -’,num2slr(i-l)]);diBp([’ ’]’) 

X 1 = inpul(*Enter the minimum x-value: *) + 1 ; 
x2 = input(’ Enter the maxumum x-valuc:*)+ 1; 

yl =input(’ Enter the minimum ’magnitude* y- value (least negative): ’) + 1 ; 

y2 = input(’Enterthe maxumum ’magnitude’ y-value (most negative): *) + l; 

yl =abs(yl);y2 = abs(y2); 

y = l:l:i;y =fliplr(y); 

yll=y(y2); 

y22=y(yl); 

clg 

cc=c(yl I:y22,xl:x2); 
mesh(cc) 

tille(’Zoomed Multiphase Energy Display of the Approximation’); 

pause 

strplt 

clg 

N = 10; 

Nl=10; 
while NK999 
N = N1; 

contour(cc,Nl,xl :x2,yl-l:y2-l) 

title(*Zoomed Multiphase Contour Display of the Approximation’); 

ylabel(’-n resolution level’) 

xlabel([’ Number of contours: ’,num2str<Nl)]) 

pause 

strplt 

Nl = input(’ Enter the number of contour levels ( < RET> to continue) ’); 
end 

ql =inpul(’Zoomin or out Further (Y/N)? ’,’s’); 

end 

clc 

% Now the Detail Signal 
dd =d(l :-lowc8l,l:LE,); 
mcsh(dd) 

tille(’ Multiphase Energy Display of the Detail’) 

pause 

strplt 

clg 

contour(dd, 10,x,yd) 

title(* Multiphase Contour Display of the Detail’); 
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ylabel(’-n resolution level’); 

xlabel(’lO contour levels’) 

pause 

strplt 

clg 

clc 

ql =input(’Doyou desire to zoom in/oul on the display (Y/N)?’,’s’); 
while ql = = ’Y’ j ql == ’y’, 

q2 = input(’Doyou want to sec the plots again (Y/N)?’.’s’); 
if q2 = = ’Y’ i q2 = = y, 

meah(dd);title(’ Previous Multiphase Energy Display of the Detail’); 
pause 

contour(dd,N,x,yd);tiUe(’PreviousMultiphase Contour Display of the Detail’); 

ylabel(’-n resolution level’); 

xlabel([ ’Number of contours: ’,num28tr(N)]) 

pause 

end 

disp([’x range: 0 to ’,num2str(j-l)]) 
disp([’y range: -1 to •’,num2str(i-l)]);dbp([ ]’) 

xl =input(’ Enter the minimum x-value:’) + 1; 
x2 = input(’ Enter the maxumum x-value:’) + l ; 

yl =input( ’Enter the minimum ’magnitude’ y-value (least negative):’); 
if yl = = 0, yl = 1; end; 

y2 = input(’ Enter the maxumum ’magnitude' y-value (most negative):’); 

yl=abs(yl);y2 = abs(y 2) ; 

y = 1 : 1 : i- 1 ;y = fliplr(y) ; 

yll=y(y2); 

y22=y(yl); 

clg 

dd=d(yll:y22,xl:x2); 

mesh(dd) 

titlc(’Zoomed Multiphase Energy Display of the Detail’); 

pause 

sUplt 

clg 

while N < 999 
contour(dd , N , x 1 : x2 ,y 1 : y2) 

litle(’Zoomed Multiphase Contour Display of the Detail’); 

ylabel(’-n resolution level’) 

xlabel([’ Number of contours: ’,num2str(N)]) 

pause 

sUplt 

N = input(’ Enter the number of contour levels (< RET > to end): ’); 
end 

ql =input(’Zoomin or out further (Y/N)? ’,’s’); 
end 

% 

% Multiple Phase Reconstruction 

mprecon 



% The End 
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%••*••••••• 

% enrgldmp.m Determining the Energy in each resolution level 
% Multi-Phase One Dimensional Case 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspiigiece. nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

% 

% For the Multi-phase case, the energy is normalized to resolution 
% level 0. All the energies of the "c’s" and the ’d’s’ of a 
% particular level should add up to the energy of the *c"s' of the 
% next higher level. The factor of two is accounted for in the Multi- 
% Phase case 
% 

% enrge — energy array of the "c’s" 

% enrgd — energy array of the ’d’s’ 

% norm - energy of the data ’cO’ 

% esum — energy sum array (enrge -1- enrgd) 

% 

% The calling routine is ’mp.m’ . 

% wavld.m and strplLm arc also necessary m-files 

% 

enrge = 2 eros(l ,-lowesl); 
enrgd = 2 eros(l , -lowest); 
a = 0; 

norm = sum(c0 .'“2); 
for row = - lowest; -1:1 
a=a-l- 1; 

enrgc(row) = sum(coeffc(row,:). *2)*2'^(-a); 
enrgd (row) = sum(cocffd( row , : ) . ^2) * 2“^ (-a) ; 
end 

enrge = enrge/nonn; 
enrge = [enrge 1]; 
enrgd = enrgd /norm; 
clg 

axis([lowesl-2 1 0 1.2]); 
subplot(21 l),bar([ lowest- 0], [enrge]) 

litJe([ ’Normalized energy of the ’c’ coefficients Multi-Phase Case’]) 

xlabelC Resolution level’) 

axis([lowest-2 1 0 1.2]); 

subplot(21 2), bar([lo west -1], [enrgd]) 

titJe( ’Normalized energy of the *d’ coefficients’) 

xlabel( ’Resolution level’) 

pause 

slrplt 

axis 



% Energy check of the previous plot 



esum = enrge + [enrgd 0] ; 
clc 

dd =[’ Verificalion of the energy in each resolution level 



Level Energy in Energy m C Difference’ 

level m-1 level m 

c+d 



]; 

disp(dd) 
j = -lowest; 
for i = 1 : -lowest 

a = csum(i) ;b = enrg c( 1 , i -I- 1 ) ; c = 1 -j;d = abs(a- b) ; 
fprintf(’%2.0f %17.15f %17.15f ’,c,a.b) 
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fprintf(’%l7.15f\n\d) 

end 

pause 

clc 

clg 

axis; 



^•***************************«****************«******* 

% mpplane.m develop the phase plane diagram 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School. Monterey CA 

% e-mail; legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

^«**************************************************** 
% Phase Plane Diagram for the One Dimensional Single Phase case 
% This routine is called by the parent routine "wavld.m*. All the 
% variables defmed in "wavld.m* are used here. New variables: 

% 

% c - ( 1 -lowest) X LL ’c* coefficient matrix 
% d - (-lowest) X LL ’d’ coefficient matrix 
% N,N1 " variables for the number of contour levels 
% 

% The main subroutine used is ptzero.m funtion m-file 

% 

clc 

clg 

% Set up the c and d matrices from the *c** and ’d** variables 

c = zero8(-lowest+ 1 , LL) ; 
d = zero8(-lowest, LL) ; 
c(l,l:length(c0)) = c0; 
clear cwork dwork 
for lvl = l:-lowest 

cw ork = cval([ ’ptzero(rmpc ’ , num2str( Ivl) . ‘ , 2 "( Ivl) - 1 ) ’]) ; 
dwork = eval([ ’ptzero( rmpd ’ , num2str( Ivl) , ’ ,2 “^ ( Ivl)- 1 ) ’ J ) ; 
c(lvl-l- 1 , 1 :length(cwork)) = cwork; 
d(lvl, 1 :length(d work)) = dwork; 
end 

c = flipud(c.*2); 
d = fUpud(d."2); 



[i j] =si 2 e<c); 

x = 0:j-l;yc=0:i-l;yd = l:i-l; 
mesh(c); 

title((’*MP* Energy Display of the Approximation Phase: ’,... 
selstr( fliplr(phB vet + 4 8))] ) ; 

pause 
N = 10; 

contour(c»10,x,yc); 

title([’*MP* Contour Display of the Approximation Phase; *.... 

setstr( fliplr(phs vet +48))]); 
ylabel(’-n resolution level’); 
xlabel(*l0 contour levels’) 
pause 
clg 
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ql =input(’Doyou dcairc lo zoom in on the display (Y/N)?\’»’); 
while ql = = ’Y' | ql = = ’y’, 
q2 = input(’Doyou want to sec the plots again (Y/N)?’,*s’); 
ifq2=-’Y* I q2 == *y’. 
mcsh(c); 

title(**MP* Previous Energy Display of the Approximation*); 
pause 

contour(c,N,x,yc); 

title(’*MP* Previous Contour Display of the Approximation’); 
ylabelC-n resolution level’); 
xlabel([ ’Number of countours: ’,num28ir(N)]) 
pause 
end 
clc 

disp([’ ’]’) 

disp([’x range: 0 to ’,num2str(j-l)]) 
disp([’y range: 0 to -’,num2str(i-l)]);disp([ ]’) 
xl =input(’ Enter the minimum x-valuc;’)+ 1; 
x2 = input(’ Enter the maxumum x-valuc:’) + l; 

yl =input(’ Enter the minimum "magnitude” y-value (least negative): ’) -t- 1; 

y2 = input(’Enterthc maxumum "magnitude" y-value (most negative);’) + 1; 

y 1 = abs(y 1 ) ;y2 = abs(y 2) ; 

y = l;l:i;y = niplr(y); 

yll=y(y2); 

y22=y(yl); 

clg 

meah(c(y 1 1 :y22 , x 1 : x2)) 

litle([’*MP* Zoomed Energy Display of the Approximation Phase: *,... 
setstr(fliplr(phsvcl +48))]) ; 

pause 

strpll 

clg 

Nl = 10; 
while NK999 
N = N1; 

contouKc(yll:y22,xl:x2),Nl,xl:x2.yl-l:y2-l) 
title([’*MP* Zoomed Contour Display of the Approximation Phase: 
setstr(flipLr(phsvct + 4 8))] ) ; 
ylabel(’-n resolution level’) 
xlabel([’ Number of contours: ’,num2str(N 1)]) 
pause 
strplt 

Nl =input(’ Enter the number of contour levels (999 to continue) ’); 
end 

ql =input(’Zoomin or out Further (Y/N)? ’.’s’); 
end 

clc 

% Now the Detail Signal 
mesh(d); 

litle([’*MP* Energy Display of the Detail Phase: 
setstr(fliplr(ph8v ct + 4 8))] ) ; 
pause 
clg 

contour(d , 1 0, x,yd) ; 

Utle([’*MP* Contour Display of the Detail Phase: ’ . 

setstr(fliplr(phsv ct + 4 8))]) ; 
ylabel(’-n resolution level’); 
xlabel(’10 contour levels’) 
pause 
clg 

clc 

ql =input(’Doyou desire to zoom in on the display (Y/N)?’, ’s’); 
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while ql = =’Y* | ql = = ’y’, 

q2 = input(’Doyou want to see the plots again (Y/N)?',’s‘); 
ifq2 = = ’Y* I q2 = = ’y\ 

mcsh(d);tilJe(’ Previous Energy Display of the Detail’); 
pause 

conlour(d.N,x.yd);tiUc(’PrcviousConlour Display of the Detail’); 
ylabel(’-n resolution level’); 
xiabel(I’ Number of contours: ’,num2slr(N)]) 
pause 
end 

disp([’x range; 0 to ’,num2sir(j-l))) 
disp([’y range: -1 to -’,num28lr(i-l)]);disp([ ]’) 

xl =input(’ Enter the minimum x-valuc:’) + 1; 
x2 = input(’Ehlerthc maxumum x-valuc:’) + 1; 

yl =input(’ Enter the minimum "magnitude" y* value (least negative):’); 
if yl = = 0, yl = 1; end; 

y2 = input(’ Enter the maxumum "magnitude" y* value (most negative):*); 

yl =abs(yl);y2 = abs(y2); 

y = l:l;i*l;y = fliplr(y); 

yll=y(y2); 

y22=y(yl); 

clg 

mcsh(d(y 1 1 :y22 . x 1 : x2)) 

title([’*MP* Zoomed Energy Display of the Detail Phase: 

8ctslr( fliplr(phs vet + 4 8))] ) ; 

pause 

strplt 

clg 

while N < 999 

contou r(d (y 11 : y 2 2 , x 1 : x2) , N , X 1 : x2 , y 1 : y 2) 
tille([’*MP* Zoomed Contour Display of the Detail Phase: 

8Ctstr(flip lr(phsvcl + 4 8))]) ; 
ylabel(’-n resolution level’) 
xlabel([’Numberof contours: ’,nura2slr(N)]) 
pause 
strplt 

N = input( ’Enter the number of contour levels (999 to end): ’); 
end 

ql = input(’Zoomin or out further (Y/N)? ’.’s’); 
cod 



function x=zOpad(Data,h,reslvl) 

^**************************** * 

% zOpad.m *** For the Multi-phase decomposition routine •** 
% One Dimensional Case 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ecc.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

% * 

% 

% This routine makes a row vector of the lower coefficients of 
% a particular resolution level. Use for only the Multi-phase 
% decomposition One-Dimensional case !! 

% 

% Data is the data vector 

% h is cither the h coeff or the g coeff vector 

% reslvl is the resolution level 

% 
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m = rcalvl; 

L=lcngth(h); 
j=L; 
x=0; 
for i = 1:L 

X = x + hO")*Ucros( 1 ,(i- l)*2"(m- l))Data zeros(l ,(j- 1 )*2"(m- 1))]; 
end 



function binary =num2bin(number, length 1) 

% 

% NUM2BIN.M Number to Binary conversion 

% 

% By: LT J. E. Legaspi 

% Profeaaor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3W4/2772 

% 

% 

% This routine conducts a conversion of a base 10 integer to a base 2 
% data array i.e. [1 0 1 0 I 1] 

% 

% The following variables are further defined; 

% 

% binary -- binary output data array 

% length 1 - input desired length of the display 

% number — base 10 integer 

% 

% The following m-files are necessary for proper operation: 

% 

% wavld.m — main program (great grandparent routine) 

% mp.m " grandparent routine 

% mprecon.m-* multiphase rcconstuction calling routine 

% 



binary = []; 
while number > .5 
number = number/2 ; 
if number- fix(number) = = .5 
number = number- . 5 ; 
binary =[1 binary]; 
else 

binary = [0 binary] ; 
end 
end 

if length(binary) < length 1 

binary = [2ero8(Labs(lengthl -length(binary)))binar>'] ; 
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function num = bin2num(bin) 

% 

% BIN2NUM.M Binary to number conversion 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam@cce.nps.navy.mil 

% Phone: (408) 646-3044/2772 

% 

% 

% This routine conducts a conversion of a base 2 data array 
% i.e. [1 0 1 0 1 1] to a base 10 integer 
% 

% The following variables are further defined: 

% 

% bin — binary output data array 
% num — base 10 integer 
% 

% The following m-filcs are necessary for proper operation: 

% 

% wavld.m — main program (great grandparent routine) 
% mp.m — grandparent routine 
% mprccon.m-- multiphase reconstuction calling routine 
% 



num = 0; 
a = length(bin); 
power = a; 
for ii = l:a 

num =num+bin(U)*2'^(powcr-ii); 
end 
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B. TWO-DIMENSIONAL ROUTINES 



^ ♦♦*•♦•***♦♦♦*♦••*♦♦•*♦•♦•♦♦ *•♦♦*♦*******♦♦♦♦*♦*•*****♦**; 

% WAV2D.M 15 SEP 92 

% Wavelet decomposition algorithm for the two-dimensional case 
% for either the single selectable (of 4 possible) or the 
% multiple phase case. 

% 



% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: leg asp i@ecc. nps.navy.mil 

% lam@ecc.nps.navy.mil 

% Phone: (408) 646-3(m/2772 

^*****»***««««***«*«*** 



% This routine does a two-dimensional dicrete wavelet decomposition 
% of a two-dimensional data array. The following routines are also 
% necessary for the proper operation of the algorithm: 



% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 



chkputer.m — checks for the minimum hardware 
mp2d.m -- 2-D multiple phase decomposition 
ptzro2d.m — 2-D row & col zero padding 
strplt.m — metafiles of selected plots 
hcoefm - routine for selected h coefficients 
phas ?? .m — routines for phase 00,01,10,1 1 decomp 
enrg2d.m — energy check for the 2-D case 
recon2d.m — 2-D reconstruction 



Variables: 



pltcnt 
Imagestr 
cO 

cl to cN 
dll todlN 
d21 to d2N 
d31 tod3N 
Hh,Hv,Gh,Gv- 
HhHv,HhGv, 

GhHv,GhGv - 2-D decomposition masks 

lowest - lowest resolution level (a neg number) 

phsvct? - phase vector for either the h or v direction 



plot counter for hard copy displays 
string variable of the input data 
resolution level 0 data input 
corresponds to the c coeff level -1 to -N 

- corresponds to the dl coefTs 

- corresponds to the d2 coeff s 
■ corresponds to the d3 coeff s 

h and g coeff s in the horiz & vert directions 



■ mm mmmmm 



% 

% initialization of the counters and plots 

clg 

axis( ’normal’); 
hold off; 
pltcnt = 0; 
chkputer 
clc 



% 

% 2-D data input routine 

Imagestr=input(’ Enter the name of the image data; ’.’s’); 
cO = eval(lmagestr) ; 

[il jl]=size(c0); 

if jl<=l I il<=l 

dispC’Bad Input Data, Restart the Program’) 
return 

end 
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clc 

% 

% 

% 2-D wavelet coefficient input routine 

hcoefin 

clg 

% 

% 

% Branch for the 2-D multiphase case 

disp((’ ’]*) 

q = input(’Doyou desire the 2-D multiphase decomposition only? ’.’s’); 
if q==y I q= = *Y’ 
mp2d 
return 
end 

% 

% 

% Initialization for the decomposition routine 
% of the single selectable phase 
% 

% initializing Cout for the decomposition routine. Cout along 
% with Dlout, D2out and D3out are intermediate working 
% variables for the determining the coefficients for a 
% particular resolution level and phase. *phas??.m" are 
% one of for possible phase decomposition routines for the 
% two dimensional case. 

Cout=cO; 

subplol(21 l),me8h(cO),tiUe(’3-DPlot of the Data Array’) 

subplol(212),conlour(cO),tJlJe(’ContouiDisplay of the Data Array’) 

pause 

strplt 

clg 

clc 

d»sp(I’ ’]’) 

q = I’Entcr the number of resolution levels desired’ 

’for decomposition of the image: ’]; 

disp(q) 

lowest = inpul( ’ : > ’) ; 
clc 

lowest = -abs( lowest) ; 

phsvcth = D; 
phsvctv = Q; 

% The actual decomposition for the single 
% selectable phase case 

for lvl = -l:-l:lowest 

clc 

q = {’Enter the desired phase decomposition: ’ 

’A) Horizontal Phase 0, Vertical Phase 0’ 

’B) Horizontal Phase 0. Vertical Phase 1 ’ 

’C) Horizontal Phase 1, Vertical Phase 0’ 

’D) Horizontal Phase 1, Vertical Phase 1 ’ 

’]; 

disp({’ ’]’) 
disp(q) 

phase = input( ’ : ’,’s’); 

if phase = = ’A’ | phase = = ’a’ 
phsvcth = [0 phsvclh]; 
phsvctv = [0 phs vetv] ; 
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phasOO 

cbcif phase = = ’B’ | phase = = ’b’ 
phsvclh = [0 phsvcth] ; 
phsvctv = [l phsvctv]; 
phasOl 

elseif phase = = ’C’ | phase = = ’c’ 
phsvcth =[1 phsvcth]; 
phsvctv = [0 phsvctv] ; 
phaslO 

elseif phase = = ’D’ | phase = =’d* 
phsvcth = [1 phsvcth]; 
phsvctv = [1 phsvctv]; 
phasl 1 

else 

disp(’ Error’) 
return 

end 

end 

% clean up for better memory utilization 

% This is primarily for 386 MATLAB 

clear wrkspe Dlout D2out D3out Gout wrk a b q qq qqq ql 

pack 

% 

% 

% Energy check for the 2-D single selectable phase case 
enrg2d 

% clean up for better memory utilization 
% This is prirnarily for 386 MATLAB 
clear E k Ivl a b c j 
pack 

% 

% 

% Display of the phase sequence in reverse bit order as per the notes 
cic 

disp((’ ’]’) 

dispC’Order of the Selected Phases from the highest to the lowest level’) 
dispG ’) 

disp( ’Horizontal: Reading left to right corresponds to going from the’) 
disp(’ higher resolution to the lower resolution level *) 
disp([’ ’]’) 
d isp(fliplr(ph8vcth)) 
disp([’ ’]’) 

disp(’ Vertical: Reading left to right corresponds to going from the’) 

dispG higherer resolution to lower resolution level ’) 

disp([’ ’]’) 

d isp ( nip lr(phsvctv)) 

disp([’ ’]’) 

disp(’<RET> to continue’) 

pause 

clc 

% 

% 

% Conduct the recomposition routine 
recon2d 

% 

clc 

disp([’ ’ ]’) 

disp(’ ROUTINE COMPLETE’) 

% The End ofWAV2D.M 
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% HCOEFIN.m 2-D wavelet coefficient input 



% 



% 

% 

% 



% By; 



% 



LT J. E. Lcgaspi 
Professor Lam 

US Naval Postgraduate School. Monterey CA 
e-mail: legaspi@ece.nps.navy.mil 



lam@ecc . nps . navy . mil 



% Phone: (408) 646-3044/2772 

% * 

% This routine obtains the desired wavelet h coefficients for both 
% the horizontal and vertical directions. The wavelets are assumed 
% separable. 

% 

% Variables: 

% 

% Hh.Hv.Gh.Gv- h and g cocffs in the horiz & vert directions 
% HhHv.HhGv, 

% GhHv.GhGv - 2-D decomposition masks 

% wavelet? • string variable in the h or v directions 

% 

% 

% Obtaining the horizontal h coefficients 
disp([’ ’)’) 

ql =[ ’Enter the number for the desired choice ’ 

’for the wavelet in the HORIZONTAL direction:’ 

’(1) Haar h coefficients ’ 

’(2) Daubechies h coefficients ’ 

’(3) Own set of h coefficients ’ 



disp(ql) 
qqq = input(’ ’); 
if qqq==3, 

Hh = inpul(’ Enter your own set of H coefficients’); 
a = sum(Hh) ;b = Hh*Hh ’ ; 
while a~=l & b~=.5 

disp(’Yourh coefficients arc NOT compactly supported!’) 

Hh = input( ’Re-enter the h coefficients or ctrl-C to stop’); 
end 

waveleth = [ ’Own’] ; 
clscif qqq==2, 

qq = inpul( ’How many taps (2-10)7 ’); 

Hh = d aubd ata(qq) ; 

waveleth = [ ’ Daubechies ’ ,num2str(qq) . ’ -tap’] ; 
cUc Hh = [.5.5]; 

waveleth = [’Haar’] ; 
qqq = l; 

end 

clc 

% 

% Now for the vertical case 
disp([’ ’]’) 

q = input(’Doyou desire the same wavelet in the vertical direction?’, ’s’); 
if q= = ’Y’ 1 q = = 'y’ 

Hv=Hh; 

wavelctv = waveleth; 
else 

disp([’ *]’) 

ql = [’Enter the number for the desired choice 
’for the wavelet in the VERTICAL direction: ’ 

’(1) Haar h coefficients 

’(2) Daubechies h coefficients ’ 

’(3) Own set of h coefficients ’ 



disp(ql) 
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qqq = input(’ '); 

if qqq = =3, 

Hv = input(’Enleryour own set of H coefficients’); 

a=sum(Hh);b=Hh*Hh’; 

while a— =1 & b~=.5 

disp(’Your h coefficients are NOT compactly supported!’) 
Hh = input(’ Re- enter the h coefficients or ctrl-C to slop’); 
end 

waveletv = [ ’Own’] ; 
clscif qqq==2, 

qq = input(’ How many taps (2-10)? ’); 

Hv =daubdata(qq); 

waveletv = [ ’ Daubechies * , num2str(qq) , ’-tap ’] ; 
else Hv = [.5.5j; 

waveletv = [’Haar’]; 

qqq=i; 

end 

clc 

end 

% 

% 

% Generating the g coefficients 

Gh = fliplr(Hh); 

for i = 2;2:length(Hh) 

Gh(l.i)=-Gh(l.i); 

end 

Gv=fliplr(Hv); 
for i=2:2;length(Hv) 

Gv(l.i)=-Gv(1.0; 

end 

% 

% 

% Generating the 2-D masks for the decomposition 
Hv=Hv*; 

Gv = Gv’; 

%the 2* factor here reduces the number of mulliplication2 for 

%each lower level coefficients 

HhHv=2*Hv*Hh; 

HhGv=2*Gv*Hh; 

GhHv=2*Hv*Gh; 

GhGv = 2*Gv*Gh; 

NHh = length(Hh); 

NHv = lengthen v); 

% End 





% PHASOO.M 2-D phase 00 decomposition of the data array 

% 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

% 

% This routine does a phase 00 decomposition given the input array - Coul. 
% horiz: 0 vert: 0 
% 

% Variables: 

% Cout, Dlout-D3out - working variables for the coefficients 

% wrkspc - woiicing array for the c coeff s 



96 



% 

% Required m-filea: 

% wav2d.m - calling program 

% decomp ll.m - plotting routine for the rea level 

% 

% 

clc 

[row8 cols]=8ize(Coul); 

% check if the row is odd, if it is, pad with one zero 
if rem(row8/2,floor(row8/2)) = = 0 
Cout=(Cout;zero8(l ,cob)]; 
row8 = rowB + 1; 
end 

%do the same with the columns 
if rcm(cols/2,floor(col8/2)) = = 0 
Cout = [Coul ze ros(rows , 1 )] ; 
cols = cols4- 1 ; 
end 

% generate the workspace 

wrkspc=[zero8(NHv-l,cols+2*NHh-3);zcro8(row8,NHh-lpout zeros(row8,NHh-2);... 

zcro8(NHv-l,cols + 2*NHh-3)]; 
clear Coul Dloul D2oul D3out 
% Now operate on the wrkspc with the decomposition masks 
a = 0; 

for rowl =1:2: rows + NHv-2 
a = a+ 1; 
b=0; 

for coll = l:2:cols + NHh-2 
b=b+l; 

wrk = wrkspc(row 1: row 1 + NHv-l , coll :coll + NHh-l); 

Cout(a,b)= 8um(8um(HhHv .* wrk)); 

Dlout(a,b)=8um(8uro(HhGv.* wrk)); 

D2out(a,b)=8um(8um(GhHv.* wrk)); 

D3oul(a,b) =8um(8um(GhGv.* wrk)); 
end 
end 

% Store the working variables with the appropriate name 

eval([’c’,num28tr(ab8(lvl)),’ = Coul;’]) 

eval([ ’d 1 * ,num28lr(ab8(lvl)) , ’ = D 1 out; *J) 

cval([ ’d2 ’ ,num28lr(ab8(lvl)) . ’ = D2oul; ’]) 

eval([ ‘d3 ’ ,num28tr(ab8(lvl)) , ' = D3out; ’]) 

clg 

% Display the variables C, Dl, D2 and D3 for the resolution level 
decomplt 



^*****************»*»*»*»***»*******************************. 

% PHASIO.M 2-D phase 01 decomposition of the data array 

% 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail; legaspi(^ece. nps.navy.mil 

% lam@ece.np8.navy.mil 

% Phone: (408) 646-3044/2772 

^***************************«*«**********************»»*****4 

% This routine does a phase 10 decomposition given the input array - Cout. 
% horiz: 1 vert: 0 
% 

% Variables: 

% Coul, Dloul- D3 out - working variables for the coefficients 

% wrkspc - working array for the c coeffs 
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% 

% Required m-fUcs: 

% wav2d.m * calling program 

% decomplt.m - plotting routine for the rea level 

% 

[rows cols] =sizc(Coul); 

% check if the row is even, if it is, pad with one zero 
if rem(rows/2,fIoor(rows/2)) = = 0 
Coul=[Cout; 2 cros(l,cols)]; 
rows = row8 + l; 
end 

%odd column check 
if rcm(cols/2,floor(cola/2)) ~ = 0 
Cout = [Cout zc ros(ro ws , 1 )] ; 
co!s = cols + l; 
end 

% generate the workspace 

wrkspc = [zeros(NHv- 1 ,cob + 2*NHh-4) ;zcros(rows.NHh-2pout zeros(row8 , NHh*2); . . . 

zcros(NHv-2,cols+2*NHh-4)]; 
clear Cout Dlout D2out D3out 
a = 0; 

for rowl = l:2:rows+NHv-2 
a = a+ 1; 
b = 0; 

for coll = l:2:cols + NHh-2 
b = b+l; 

wrk = wrkspc( row 1 : row 1 + NK v- 1 , col 1 : col 1 + NHh- 1 ) ; 

Cout(a,b)= 8um{sum(HhHv .* wrk)); 

Dlout(a,b) = 8um(8um(HhGv.* wrk)); 

D2out(a,b)=8um(8um(GhHv.* wrk)); 

D3out(a,b) = 8um(8um(GhGv.* wrk)); 
end 
end 

lvl2 = abs(lvl); 

cval([’c’,num2str(lvl2),* = Cout;*]) 
cval([ *d 1 ’ ,num2str(lvl2) , ’ = D 1 out; ’]) 
cval([’d2’ ,num28tr(lvl2) , ’ = D2out; *]) 
cval([’d3 ’ ,num2str(lvl2) , ’ = D3oul; *]) 
clg 

decomplt 



% 

% PHASOl.M 2-D phase 01 decomposition of the data array 

% 

% 

% By: LT J. E. Lcgaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: lcgaspi@ecc.np8.navy.mil 

% lam@ccc.nps.navy.mil 

% Phone: (408) 646-3044/2772 

%**•**••*•*•****♦ * * 

% This routine docs a phase 01 decomposition given the input array - Cout. 
% horiz: 0 vert: 1 
% 

% Variables: 

% Cout, Dlout-D3out - working variables for the coefficients 

% wrkspc - working array for the c coeff s 

% 

% Required m- files: 

% wav2d.m - catling program 

% decomplt.m - plotting routine for the res level 
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% 

clc 

[rows cols] =8izc(Coul); 

% check if the row is odd, if it is, pad with one zero 
if rem(row8/2,floor(row8/2)) - = 0 
Coul= [Coul;zcro8( I .cols)] ; 
rows = rows + 1 ; 
end 

% check if the coluinns arc even 
if rcm{col8/2,floor(cols/2)) = = 0 
Cout=[Cout zcro8(row8,l)]; 
cols = cols + 1; 
end 

% generate the workspace 

wrkspc = [zcro8(NHv-2,cols+2*NHh-3);zero8(row8.NHh-l)Coul zcro8(row8,NHh-2);. .. 

zcro8(NHv-2,cols + 2*NHh-3)]; 
clear Coul Dlout D2out D3out 
a = 0; 

for rowl = 1:2: rows + NHv*2 
a = a+l; 
b = 0; 

for coll = l;2;cob + NHh-2 
b=b+l; 

wrk = wrk8pc(rowl:rowl + NHv-l, coll: coll + NHh-l); 

Cout(a,b)= 8um(8um(HhHv . • wrk)); 

Dloui(a,b)=8um(8um(HhGv.* wrk)); 

D2out<a,b)=8um(8um(GhHv.* wrk)); 

D3oul(a,b)=8um(8uin(GhGv.* wrk)); 
end 
end 

lvl2 = abs(lvl); 

cval([ ’c’ ,num28lr<lvl2) , ’ = Cout; ’]) 
cval([’d r,num28tr(lvl2) ,’ = Dlout; ’]) 
cval([’d2’,num28tr(lvl2).’ =D2oul;’]) 
cval([’d3’,num28lr<lvl2),’ =D3out;’]) 
cig 

decompit 



% ******************* 

% PHASl l.M 2-D phase 1 1 decomposition of the data array 



% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: lcgaspi@ece.np8.navy.mil 

% lam@cce.nps.navy.mil 

% Phone: (408) 646-3044/2772 

% * 

% This routine does a phase 1 1 decomposition given the input array - Cout. 
% horiz: 1 vert; 1 
% 

% Variables: 

% Cout, Dlout-D3out - working variables for the coefficients 

% wrkspc - working array for the c cocfC s 

% 

% Required m-filcs: 

% wav 2d. m - calling program 

% decompit. m - plotting routine for the res level 

% 

% 

clc 
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[rows colsj =8ize(Coul); 



% check of the row is odd, if it is, pad with one zero 
if rem(rows/2,noor(row8/2)) ~ = 0 
Cout = [Coui;zcro8( 1 ,cols)] ; 
rows = rows + 1 ; 
end 

%do the same with the columns 
if rem(cob/2,floor(cols/2)) — = 0 
Cout = [Cout zcro8(row8,l)]; 
col8=col8+ 1; 
end 

% generate the workspace 

wrkBpc = [zero8(NHv-2,cols+2*NHh-4);zcro8(row8.NHh-2poui zcros(row8,NHh-2);. . . 

zcro8(NHv-2,cob+2*NHh-4)J; 
clear Cout Dlout D2out D3out 
a = 0; 

for rowl = l:2:rows+NHv-2 
a=a+ 1; 
b=0; 

for coll = l:2:cob + NHh-2 
b = b + l; 

wrk=wrks|>c(rowl;rowl +NHv*l,coll:coll +NHh-l); 

Coui<a,b)= 8um(sum(HhHv .• wrk)); 

Dloul(a.b)=8um(8um(HhGv.* wrk)); 

D2oul(a,b) = 8um(sum(GhHv.* wrk)); 

D3out(a,b) = 8um(8um(GhGv.* wrk)); 
end 
end 

lvl2=abs(lvl); 

cval([ ’c’ ,num28tr<lvl2) , ’ = Cout; ’J) 
cval([’dr,num2str(lvl2),’ = Dlout;’]) 
eval( [ ’d2’ ,num2sti<lvl2) , * = D2out; ’]) 
cval( [’d3 ’ ,num28tr(lvl2) , ’ = D3out; ’]) 
clg 

decomplt 





% DECOMPLT. M dbpby8 the coefficients for a particular level 

% 

% By: LT J. E. Legaspi 

% Professor Lara 

% US Naval Postgraduate School, Monterey CA 

% e-mail: lcgaspi@ccc.np8.navy.mil 

% lam@ecc.nps.navy.mil 

% Phone: (408) 646-3044/2772 



% Variables: as defined in WAV2D.M and PHAS77.M 

% 

% Required m-files: 

% WAV2D.M parent routine 

% PHAS77.M one of four possible calling routines 

% STRPLT.M metafile storage of selected plots 

% 

clg 

subploi(22 1 ) , conlour(Cout) ; xbbcl( ’ Cm’ ) ; 
tiUc([ ’wavelet horiz: ’.wavclcth]); 

8ubploi(222) , conlour< D lout); xbbcl( ’Dim’) 
litJe([ ’wavelet vert: '.wavclctvj) 
subplol(223 ) , contour(D2out) ; tiUc( ’ D2m ’) 
xiabcl([’ Horiz phase: \8Ctstr(niplr(phsvcth-t-48)).’ Vert phase: 
8ctstr(fliplr(phsvctv +48))]) 
subploi(224),contour(D3out);tiUc(’D3m’) 
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xlabel([’Resolulion Level ’,num28lr(lvl)]) 

pause 

slrpll 

clg 

Bubp Iol(22 1 ) , mc8h(Cout) ; xlabel{ ’Cm ’) 
liUe{[’wavclet horiz: ’,wavelc!h]); 
subplot(222) , mc8h(D 1 out) ; xlabel{* D Ira’) 
uUe([ ’wavelet vert: ’.waveletv]) 
gubplot(223) , mesh(D2out) ; tiUe( ’ D2m ’) 

xlabel([ ’Horiz phase: ’,8eUtr(niplr(phsvcth+48)),’ Vert phase; 

8elstr(niplr(phsvctv +48))]) 

8ubplot(224) ,me8h(D3out);tiUc( ’D3m’) 
xlabel{[’ Resolution Level ’,num28tr(lvl)]) 
pause 
strplt 

% END 



% 

% enrg2d.m Determining the Energy in each resolution level 
% Two Dimensional Case 

% 

% By: LT J. E. Lcgaspi 

% Professor Lara 

% US Naval Postgraduate School, Monterey C A 

% e-mail: legaspi@ece.np8.navy.mil 

% lam@ece.np8.navy.mil 

% Phone; (408) 646-3CW4/2772 



% For the single phase case, the energy is normalized to resolution 
% level 0. All the energies of the "c’s* and the "d’s" of a 
% particular level should add up to the energy of the ‘’c’s* of the 
% next higher level. The calling routine is "wavTd.m* and ’’strplt.m" 

% is a necessary m-file. 

% 

% Variables (other than defined in mp2d.ra): 

% 

% E - energy array for [c,dl,d2,d3] for caach level 

% csum - total energy of each level: [c(m) + dl +d2+d3]=c(m+ 1) 

% 

% Initialization of the energy sum for resolution level 0 
E=zeros(- lowest + 1,4); 

E( 1 . 1 ) = 8um(8um(c0 . "2)) ; 

lvl=-l; 

k=l; 

% Calculate the energies for each resolution level 
while lvl> = lowest 

E(k+ 1 ,1) =8um(8um((eval([’c’,num28Ur(ab8(lvl))])> *2)); 
E(k+l,2)=8um(8um((eval([’dl ’,num2slr(abs(lvl))]))'’‘2)); 

E(k+ 1 ,3) =8um(8um((eval([’d2’,num28lr(ab8(lvl))])> *2)); 

E(k+ 1 ,4) = 8um(8um((eval([’d3 ’,num28tr(ab8(lvl))])> "2)); 
lvl = lvl-l; 
k = k+l; 
end 

E=E/E(1,1); % normalize the energy to level 0 

E=flipud(E); % top row is the lowest resolution 
% Display of the energies by coefficients 
clg 

k = [lowesM 1 0 1.1]; 
axis(k);subplot(221),bar(lowcst:l:0,E(:,l)) 
xlabel(’c’) 

liUe([ ’wavelet horiz: ’.waveleth]) 



101 



axis(k);subplol(222),bar(lowc8l: 1 :0,E(: ,2)) 
xlabel(’dl ’) 

tilJe([ ‘wavelet vert; ’.waveletv]) 
axis(k);8ubpIot(223),bar(lowest: 1 :0,E(:,3)) 
tiUc(’d2’) 

xlabei( ’single phase case’) 
axis(k) ;subplot(224) ,bar( lowest: 1 : 0, E( : ,4)) 
tiUc(’d3’) 

pause 

slrplt 

cig 

clc 

E=E’; 

esum = sum(E); 
diflp([’ ’]’) 

dd=[’ Verification of the energy in each resolution level ’ 

’Level Energy in Energy in C Difference’ 

’ m level m-1 level m ’ 

c+dl+d2+d3 

’]; 

disp(dd) 
j = -lowest; 
for i = l:-lowest 

a = esum(i) ; b = E( 1 , i + 1 ) ; c = - j + 1 ; d = abs(a-b) ; 
fprintf(’%2.0f %17.15f %17.15f ’.c.a.b) 

fprintf{’ %l7.l5fVn’,d) 

end 

pause 

clc 

% End 



% recon2d.m 2-D reconstruction algorithm for the single 
% selectable phase case 

% 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

% 

% Variables (other than defined in wav2d.m): 

% 

% HvHh.HvGh - rcconstrution masks 
% GvHh.GvGh 

% whichl - phase decomposition for the particular level 

% cwork - reconstructed c coefficient array for a level 

% 

% Required routines: 

% wav2d.m 

% strplt.m 

% mmlt??.m - ?? = 00,01 ,10,1 1 for the 

% mask multiplication 

% 

% 

clc 

disp([’ ’]’) 

q = input(’Doyou desire to do the reconstruction (Y/N)? ’,’s’); 
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if q = =’y’ |q==:’Y’ 

% The following four lines generate the reconstruction masks 
H vHh = niph< Hipud (HhH V)) ; 

HvGh = niplr<flipud(GhHv)) ; 

GvHh = fUplr<fUpud(HhGv)); 

GvGh = niplr(flipud(GhGv)); 
disp([’ ’]’) 

q = inpul(’Do you desire to see the masks? ’,’s’); 
if q==‘Y’ I q==’y’ 

% We will border the arrays with zeros for a better mesh plot 
[aa bb] =size(HhHv); 
dispC ’) 

dispC’HhHv decomposition’) 
disp(HhHv) 
subplol(121) 
axis(’ square’) 

mesh([zeros( 1 , bb + 2) ; zeros(aa , 1 )HhH v zeros(aa , 1 ) ; zeros( 1 ,bb + 2)] ) 
litle(’HhHv decomposition mask’) 
disp(* ’) 

disp(’HvHh reconstruction*) 
disp(HvHh) 
subplot(122) 
axis( ’square’) 

mesh([zeros(l ,bb + 2) ;zeros(aa, l)HvHh zeros(aa, 1) ;zeros( 1 ,bb + 2)]) 
litie(*HvHh reconstruction mask’) 
pause 
strplt 
clg 

disp(’ ’) 

disp(’HhGv decomposition’) 
disp(HhGv) 
subplot(121) 
axis( ’square’) 

mesh([zeros(l,bb + 2);zero8(aa,l)HhGv zeros(aa.l);zeros(l,bb + 2)]) 
litle(’HhGv decomposition mask’) 
disp(’ ’) 

disp(’GvHh reconstruction’) 
disp(GvHh) 
subplot(122) 
axis(’square’) 

mesh([zeros(l ,bb + 2);zeros(aa, l)GvHh zeros(aa,l);zeros(l ,bb+2)]) 
tiUe(’GvHh reconstruction mask’) 
pause 
strplt 
clg 

disp(’ ’) 

disp(’GhHv decomposition’) 
disp(GhHv) 
subplot(121) 
axis( ’square’) 

mesh( [zer os( 1 ,bb + 2) ; zeros(aa, 1 ) G hH v zeros(aa , 1 ) ; zeros( 1 , bb + 2)]) 
liUef’GhHv decomposition mask’) 
disp(’ ’) 

disp(’HvGh reconstruction’) 
disp(HvGh) 
subplol(122) 
axis(’square’) 

mesh([zero8(l,bb + 2);zeros(aa,l)HvGh zeros(aa,l);zeros(l.bb + 2)]) 
litle(’HvGh reconstruction mask’) 
pause 
strplt 
clg 

disp(’ ’) 

disp(’GhGv decomposition’) 
disp(GhGv) 
subplot(121) 
axis( ’square’) 
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roe8h([zcro8{ 1 , bb + 2) ;zcro8{aa , 1 )G hG v zcro8(aa , 1 ) ; zcro8( 1 , bb 4- 2)] ) 
tiUe(’GhGv dccompoailion mask*) 
dispC ’) 

disp(’GvGh reconstruction’) 
disp(GvGh) 

8ubplot(122) 

axisCsquarc’) 

mesh([zcro8( 1 ,bb + 2) ;zero8(aa. 1 )G vGh zeroa(aa, 1 ) ;zeroa( 1 ,bb + 2)]) 
tiUc(’GvGh reconstruction mask’) 
pause 
strplt 
clg 
end 

clear aa bb 
axis(’ normal’) 
clc 

disp([’ ’]’) 

dispC NOW DOING THE RECONSTRUCTION’) 

%-- 

% start at the lowest level 
% and calculate the c’s for the next level 
lowest = -ab8(lowest) ; 
lvl = lowest; 

cwork =eval([’c* ,num28tr(ab8( lowest))]); 
aa = 0; 
while IvKO 
aa==aa + l; 

[a b] =size<eval([’c’,num28tr(ab8(lvl+ 1))])); 

dl =eval([’dl’,num28tr(ab8(lvl))]); 

d2 = eval([’d2’,num28tr<ab8(lvl))]); 

d 3 = eval( [ ’d3 ’ ,num28tr(ab8(lvl))] ) ; 

whichl = [num28tr(phsvcth(aa)) ,num28lr(phsvctv(aa))] ; 

eval([’lc =mmlt’,whichl, ’(cwork, HvHh);’J) 

eval([ ’ td 1 = mmlt’ .which 1 , ’ (d 1 ,G vHh) ; ’]) 

cval([’td2=mmlt’, whichl, *(d2,HvGh);’]) 

eval([’td3 = mmlt’, whichl , ’(d3 ,GvGh); ’]) 

clear cwork dl d2 d3 

cwork = tc(l:a.l:b)+tdl(l:a,l-.b) + td2(l:a.l:b)+td3(l:B.l;b); 

clear tc tdl td2 td3 a b whichl 

clg 

axiaCsquare’) 

% plot the result of the c’s for comparison 
8ubplot( 1 2 1 ) ,conlour(cwork) ; 
title( ’ Reconstructed C array’) 

8ubplot(122),contour(eval([’c’,num28tr(ab8(lvl)-l)])); 
title([ ’Actual C array’]) 
xlabel([’ Resolution Level *,num2str(lvl+ 1)]) 
pause 
strplt 
clg 

8ubplot( 121) .mesh(cwork) ; 

titJe(’ Reconstructed C array*) 

8ubplot( 1 22) ,mesh(eval([ ’c ’ ,num2stj<ab8(lvl)- 1)])); 
titlc([ ’Actual C array’]) 
xlabel([ ’Resolution Level ’,num2str(lvl+ 1)]) 
pause 
strplt 
clg 

errr = cwork- eval([ ’c ’ , num2str<ab8(lv 14-1))]); 
if max(max(abs(enT))) — = 0 
mesh(ab8(errr)) 

titlc([ ’Absolute error from ’,num2str(lvl),’ to ’,num2str(lvl4- 1)]) 
xlabcl([’ Maximum Error Value: ’,num2str(max(max(ab8(errr))))]) 
pause 
strplt 
else 
clc 
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Clg 

disp([’ ’]’) 

dispC ••• NO ERROR IN THE RECONSTRUCTION *•*’) 

disp([’ ’]’) 

dispC NOW CALCULATING THE NEXT RECONSTRUCTION’) 

disp([’ reOM LEVEL ’num28lr(lvl+l),’ TO 

num2slr(lvl+2)]) 

end 

lvl = lvl+l; 
end 

ax is( ’normal’) 

clc 

end 



function x = 

^ ****••*< 



mmltOO(airay,maak) 



% MMLTOO.M mask mulliplicalion and summation for the 2-D phase 00 
% reconstruction 

% 

% By; LT J. E. Legaspi 
% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: lega8pi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

^***********************»********»**************************« 

% !!!!.»!!! This routine is only for PHASE-00 !!!!!!!!! 

% Required M = files : 

% wav2d.m - grandparent routine 

% recon2d.m— parent routine 

% ptzro2d.m — puts n zeros rows & cols btwn coeff s 

% 



array = ptzro2d(airay , 1 ) ; 

% generate the workspace 
[NHv NHh]=sizc(maflk); 

[row col] =size(array); 

array = [array zeros(row,NHh-2)] ;cols = col+ NHh-2; 
array = [array; zero8(NHv-2,cols)]; 
a = 0; 

for row 1 = 1: row- 1 



a = a-t- 1 ; 

b = 0; 

for coll = l:col-l 
b = b+l; 

wrk = array (row 1 : row 1 -I- NHv- 1 ,coll : col 1 + NHh- 1 ) ; 
x(a,b) = sum(sum(mask . * wrk)) ; 
end 
end 
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function x = mmU10(arTay,mask) 

^ ***••*•*•***•••*•*••*•*•*****»**** *•*••*••**••«■***••«•*•*• ************ 

% MMLTIO.M mask multiplication and summation for the 2-D phase 10 
% reconstruction 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps.navy.mil 

% lam@cce.nps.navy.mil 

% Phone: (408) 646-3044/2772 

%**••** * * 

% !!!!!!!! This routine is only for PHASE-10 !!!!!!!!! 

% Required M= files: 

% wav2d.m -- grandparent routine 

% recon2d.m-- parent routine 

% ptzro2d.m -- puts n zeros rows & cols btwn coeff s 

% 

array = ptzro2d (array , 1 ) ; 

% generate the workspace 
(NHv NHh] =size(mask); 

[row col]=size(aiTay); 

array = (zeros(row, 1) array zeros(row,NHh-2); zero8(NHv-2,col + NHh-1)]; 
a=0; 

for rowl = l:row-l 
a = a-f-l; 
b = 0; 

for coll = l:col 
b = b + l; 

wrk = array(row 1: row 1 + NHv- 1, coll :colH- NHh-1); 
x(a,b) = 8um(sum(mask . * wrk)) ; 
end 
cod 



function x = mmltOl (array .mask) 



% MMLTOl.M mask multiplication and summation for the 2-D phase 01 
% reconstruction 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ecc.nps.navy.mil 

% lam@ecc.nps.navy.mil 

% Phone: (408) 646-3044/2772 



% !!!!!!!! This routine is only for PHASE-01 !!!!!!!!! 

% Required M = files: 

% wav 2d. m — grandparent routine 

% rccon2d.m — parent routine 

% ptzro2d.m— puts n zeros rows & cols btw'n coeffs 

% 

array = ptzro2d(airay , 1) ; 

% generate the workspace 
[NHv NHh] =sizc(maBk); 

[row col] = sizc(aiTay) ; 

array = [zeros(l ,col+NHh-2);arrayzeros(row,NHh-2);zeros(NHv-2,coH-NHh-2)]; 
a = 0; 

for rowl = l:row 
a = a+ 1 ; 
b = 0; 

for coll = 1 :col-l 
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b=b+ 1 ; 

wrk = array(rowl:rowl + NHv-l, coll.-coil +NHh-l); 
x(a,b) =8um(8um(mask. *wrk)); 
end 
end 



foncUoD x = mmlll 1 (array, mask) 

^, ********************************************************************* 

% MMLTl 1 .M mask multiplication and summation for the 2-D phase 1 1 
% reconstruction 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School. Monterey CA 

% e-mail: legaspi@ece.np8.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

^*****«******««***«*«******«««******«*«********«*********************** 

% !!!!!!!! This routine is only for PHASE-11 !!!!!!!!! 

% Required M = files: 

% wav2d.m — grandparent routine 

% recon2d.m-- parent routine 

% ptzro2d.m-- puts n zeros rows & cols btwn coeffs 

% 

array =ptzro2d(array , 1 ); 

% generate the workspace 
[NHv NHh)=size(ma8k); 

[row col] =8ize(array); 

array =[zero8(l ,col);array); rows = row + 1 ; 

array = [zcros( rows , 1 ) array] ; CO b = col + 1 ; 

array = (array zeros(row8, NHh-2)] ;cob = cob + NHh-2; 

array =(array; zeros(NHv-2,cob)]; 

a = 0; 

for rowl = l:row 
a=a+ 1; 
b=0; 

for coll = l:col 
b = b+l; 

wrk = array( row 1 : row H-NHv-l.coll: coll -1-NHh-l); 
x(a,b) = 8um(sum(mask . * wrk)) ; 
end 
end 

function x = ptzro2d(input,n) 



%PTZR02D.M puts *n" rows and columns of zeros after each 
% coefficient in the 2>D input array 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: lcgaspi@ece.nps.navy.mil 

% lam@ecc.nps.navy.mil 

% Phone: (408) 646-3044/2772 

^*»*»*********«*,*******,**************************************** 
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% Required routines: 

% wav2d.m -- great-grandparent routine 

% recon2d.m — grand parent routine 

% mmJt??.m -- parent routine, ?? =00, 01, 10,1 1 

% 

[r c] =size(input); 

xO=D; 

for j = l:c 

x0 = [x0 inpul(:,j) zero8(r,n)]; 
end 

(r cj = 8ize(x0); 

x=D: 

for j = l:r 

X = [x;x0(j, :);zcro8(n,c)] ; 
end 



^••••******************** 

% MP2D.M 2- D Multi-phase decomposition 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.nps. navy. mil 

% lam@ece.np8.navy.mil 

% Phone: (408) 646-3044/2772 

% 

% MP2D conducts a multiple phase decomposition routine for the 2-D 
% case. The following m-files are required: 

% 

% wav2d.m - calling routine 

% mpdcmp2d.m - multiple phase decomp for 1 resolution level 
% strplt.m - metafile storage of selected plots 
% enrg2dmp.m - multiple phase energy check 

% 

% Variables are the same as the calling program 

% 

% 

% 

% Initialization of the routine 

clc 

cig 

subplot(21 l),mcsh(c0),tille(’3-DPlot of the Image Array’) 

subplot(212),contour(cO),Utle<*ConlourDisplay of the Image’) 

pause 

strplt 

cIg 

cIg 

axis(’ normal’) 
hold off 
disp([’ ’]’) 

lowest =input(’ Enter the lowest resolution level; ’); 
lowest = -abs(lo west); 
m = 0; 

% 

% 

% Calculating the Coefficients for each resolution level 
while m > = lowest + 1 

cval([ ’c ’ ,num2str(-m + 1 ) , ’ = mpdcmp2d(c ’ , num2str(ab8(m)) , ’ , H hH v ,m) ; ’] ) 
eval([ ’d r ,num28tr(-m + 1 ) , ’ = mpdcmp2d(c ’ , num2str(abs(m)) , ’ , HhG v , m) ; ’)) 
cval([ ’d2’ ,num28tr(-m + 1), ’ = mpdcmp2d(c’ , num2str(abs(m)), ’ , GhH v,m); ’]) 
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eval([’d3’,num28tr^'m+ 1)/ =mpdcmp2d(c’,num28lr(abs(m))/,GhGv,m);’]) 
% Display of the coefficients 

8ubplol(221 ) ,eval([ ’mesh(c ’ ,num28tr(ab8(m) + 1 ) , ’) ’]) 
xlabel([’c*]) 

tiUe([’wavelet horiz: ’.waveleth]) 

8ubplot(222),cval([*mesh(dl \num28lr(abs(m)+ 1),’)’]) 
xlabcKI’dl’]) 

tiUe([ ’wavelet vert: ’.wavcletv]) 

8ubplol(223) ,cval([ ’mesh(d2 ’ ,num28lr(abs(m) + 1 ) . ’) ’]) 
tiUe((’d2*]) 

xlabel( ’multiphase case’) 

8ubplol(224) ,cval([’mesh(d3 ’ ,num28lr(ab8(m) + 1 ), ’) ’]) 
tiUe([’d3’]) 

xlabel([ ’Resolution level ’,num28tr(m-l)]) 
pause 
8lrplt 
clg 

subplol(22 1 ) ,eval([’conlour<c’ ,num28tr(ab8(m) + 1 ) , ’) ’]) 
xlabel([’c’]) 

tiUe([ ’wavelet horiz: *, waveleth]) 

8ubplot(222) ,cval([’contour(d 1 ’ ,nura2atr(ab8(m) + 1 ) , ’) ’]) 
xlabcl([’dl’]) 

tiUe([ ’wavelet vert; *,waveletv]) 
subplot(223) ,cval([’contour(d2 ’ ,num2str(ab8(m) + 1 ) /) ’]) 
uUe([’d2’]) 

xlabel( ’multiphase case’) 

8ubplol(224) ,eval([’contour(d3 *,num28tr(ab8(m) + 1 ). *) ’]) 
tiUe([’d3’]) 

xlabel([’ Resolution level ’,num28lr^m-l)]) 
pause 
strplt 
clg 

m=m-l; 

end 

% 

% 

% Energy check for the 2-D multiphase case 
enrg2dmp 

% 

% The End 



function X = mpdcmp2d(array .mask ,ra) 

^***«********************«*«**********************«********* •••«*«** 

% MPDCMP2D.M 2-D Multi-phase decomposition decomposition function 
% m-file 

% 

% By: LT J. E. Legaspi 

% Professor Lara 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ece.np8.navy.mil 

% lam@ece.np8.navy.mil 

% Phone: (408) 646-3044/2772 

^**********************,******************************************** 

% MPDCMP2D is the actual function m-file that zero pads the input 
% image array given the size of the mask operator and the 

% resolution level m. "m" is the current resolution level 

% and is a negative number. 

% 

% mp2d.m - calling routine 



%■ 
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[NHv NHh] =8izc(maflk); 

[a b]=8izc(array); 

x = zcro8(a + 2*(-m)*(NHv-l),b + 2"(-m)*(NHh-l)); 
fork = l:NHh 
for l = l:NHv 

%add the row8 of zcro8 

work = [zcro8(2^(-m)*(l-l),b);array;zcro8((NHv-l)*2"(-in),b)]; 
[a b] =8izc(work); 

%add the cols of zeros 

workl =[zcro8(a,2'"(-m)*(k-l))woiit zcros(a,(NHh-k)*2^(-m))]; 
work2 = mask(NHv-H-l,NHh-k+l)*workl ; 
x = x + work2; 

% dear work workl work2 
end 
end 



% cnrg2dmp.m Determining the energy in each resolution level 
% multiphase case 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: legaspi@ccc.nps.navy.mil 

% lam@ece.np8.navy.mil 

% Phone: (408) 646-3044/2772 

% In this MULTIPHASE ease, the lower level will have 
% four possible phases, thus its energy will be 4 times as high as 
% the next higher level so we will divide the energy of the lower 
% level and the value should equal the energy of the c’s of the 
% next higher level. 

% 

% Variables (other than defined in mp2d.m): 

% 

% E - energy array for [c,dl,d2,d3] for eaach level 
% esum - total energy of each level: [c(m)+dl +d2+d3J/4=c(m+l) 

% 

E =zcros(- lowest -t- 1 ,4); 

% Initialization of the energy for resolution level 0 
E( 1 , 1) = 8um(sum(c0 . *2)) ; 
lvl = -l; 
k = l; 

% Calculate the energies for each resolution level 
while lvl> = lowest 

E(k+l,l) = 8um(8um((cval([’c’,num2stj<ab8(lvl))])).*2*4*(lvl))); 

E(k -1-1,2) = 8um(8ura((cval([’d 1 ’ ,num2str(ab8(lvl))])) ''2*4'^(lvl))) ; 
E(k-fl,3) = 8um(8um((cval([’d2’,num28U<ab8(lvl))]))"2*4"(lvl))); 

E(k-H 1 ,4) = 8um(8um((cval([’d3 ’,num28tr(ab8(lvl))])) '^2*4'^(lvl))); 
lvl = lvl-l; 
k=k-l-l; 
end 

E=E/E(1 ,1); % normalize the energy to level 0 
E=flipud(E); % top row is lowest resolution level 
% Display of the energies by cocfficienls 
clg 

k = [lowcsM 1 0 1.1]; 
axis(k) ;subplot(22 1 ) , baKlowest: 1 : 0 , E( : , 1 )) 
xlabcl(’c’) 
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liUc( [ ’wavelet horiz; ’,wavelcth]) 
axifl(k);8ubplot(222),bai<lowc8f 1 :0,E(:,2)) 
xlabcK’dl’) 

tiUe([’ wavelet vert: ’.waveletv]) 
axis(k);8ubplot(223),bar(lowc«t: 1 :0.E(: ,3)) 
tiUe(*d2’) 

xlabcl(’multiphasc case’) 
axis (k);8ubplot(224),bar( lowest: 1 :0,E(: ,4)) 
tiUe(’d3’) 

pause 

slrplt 

clg 

clc 

E=E’; 

esum=sum(E); 
e8um=8um(E); 
disp([’ ’]’) 

dd = [’ Verification of the energy in each resolution level 

’Level Energy in Energy in C Difference’ 

’ m level m- 1 = level m ’ 

* c+dl+d2+d3 

‘ ]; 

disp(dd) 
j = lowest; 
for i = l:-lowest 

a - esum(i);b = E(1 ,i+ l);d =ab8(a-b);c =j + 1 ; 
fprintf(’%2.0f %17.15f %17.15f ’.c.a.b) 

fprintf(’%17.15An’,d) 

j=j+i; 

end 

pause 

clc 

% End 



% IDATA.M Sample Image data for testing the routines 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey C A 

% e-mail: legaspi@ece.np8.navy.mil 

% lam@ece.np8.navy.mil 

% Phone: (408) 646-3044/2772 

^**********«******«****«*««*******«******«*«******«**********«*****»«*«** 

K = mcnu( ’Select image’ , ’Boxl ’ , ’Box2 ’ , ’Box3 ’, ’B/W’ , ’2B/W’ , ’ 8B/W’ .... 

’X’, ’Hexagon’); 



if K= = l 

% first image is a box 
im=zcro8(100,100); 
im(20:21, 20:80) =ones(2.61) 
im(79:80.20:80) = ones(2.61) 
im(22:78,20:2^1) =oncs(57,2) 
im(22:78.79:80) = oncs(57,2) 



elscif K==2 
im = zcro8( 1 00, 1 00) ; 
im(20,20:79) = ones(1.60); 
im(80,20:79)=oncs(l,60); 
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im(20:79,20) = onc8(60,l); 
im(20:79.80)=ones(60,l); 

elscif K==3 
im = zeros( 1 00 , 1 00) ; 
im(l:2,l:61)=ones(2,61); 
im(60:61,l:61) = oncs(2,61); 
im(3:59,l:2) = ones(5 7,2); 
im(3:59, 60:61) =ones(57,2); 

elseif K==4 

% second image is a 1 white and 1 black rectangle 
im = 2 cro 8 ( 1 00, 1 00) ; 
im(:,51:100)=ones(100,50); 



ciseif K==5 

% third image arc two black boxes and two while boxes 

im = zeros( 100,1 00) ; 

im(l :50, 1 :50) =oncs(50,50); 

im(5 1 : 1 00, 5 1 : 1 00) = oncs(50, 50) ; 



cbcif K= = 6 

% fourth image arc 8 black and 8 white boxes 
im=zeros(100,100); 

im(l:25,;) =[ones(25,25)zeros(25,25) oncs(25,25) zeros(25.25)]; 
im(26:50,:) = [zcro8(25,25)ones(25,25) zeros(25,25) oncs(25,25)]; 
im(51:75,:) = [ones(25,25)zcro8(25,25) oncs(25,25) zcros(25,25)]; 
im(76:l00,:)=(zero8(25,25)ones(25,25) zcros(25,25) ones(25,25)]; 

elscif K= =7 

% fifth image are 2 45 degree diagonals 

im=eyc(l00,l00); 

im =aign(im + fUpud(im)); 



elscif K= =8 

% sixth image is a hexagon 
for row =1:20 
n=2*row-l; 

im6(row,n:n + 2) =ones(l ,3); 
end 

[a b] =size(im6); 

row2 = ceil(sqrt(a*2 + b*2)) ; 

im6a = zeros(row2 ,b) ; 

im6a(: ,b-l :b) =oncs(row2,2); 

rim = [im6 ; im6a ; fl ip M im6)] ; 

lim = fUplr(rim); 

im6 = [lim rim] ; 

clear row row2 im6a lim rim 

[a b]=size(im6); 

im6a = zcros( 1 00 , 1 00) ; 

aa = floor(( 1 00-a)/2) ; 

bb = floor((100-b)/2); 

im6a(aa + 1 : aa + a , bb + 1 : bb + b) = iro6 ; 

im = im6a; 

clear im6 im6a a aa b bb K n 
end 
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C. ROUTINES COMMON TO BOTH DIMENSIONS 



%••••• 

% chkputer.m checks for the proper ope rating system 

% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: lcgaspi@ece.nps.navy.mil 

% lam@ece.nps.navy.mil 

% Phone: (408) 646-3044/2772 

^««««*«******************«*«****«*»**********«»*»»»»,*** 

% 

% This routine reads the variable "emptr" to sec if it can run in 
% the current system. At this time, the algorithm will run on a 
% PC386, UNIX, or SUN system. If the operating system is a VAX- 
% based, the routine will not be as reliable since the VAX MATLAB 
% docs not recognize IEEE arithmetic. The routine also deletes all 
% other leg*. met files prior to running. 

% 

[emptr mxsz] =computer; 
a = lenglh(cmptr); 
clc; 

if cmptr(l) = = ’P’&cmptr(2)= =’C’&cmptr(3) = = ’3’, 

!del leg*. met; 

clseif cmptr(l) = = ’S’&cmptr(2)= =*U’&cmptr(3) = = *N’, 

!nn leg*. met; 

elseif cmptr(a-3)= = *U’&cmpU<a)= =’X*. 

!rm leg*. met; 

else 

clc 

disp([’ ’]’) 

ql =['*****Errorin input******’ 

’Hardware cannot support ’ 

’the algorithm... PLS ’ 

’ Restart the Program *]; 
disp(ql) 
return 
end 



% strplt.m 
% 

% Store Plot M-filc for wavelet decomposiUon routine 
% 

% By: LT J. E. Legaspi 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% e-mail: Iega8pi@ece.nps.navy.mil 

% lam@ecc.nps.navy.mil 

% Phone: (408) 646-3044/2772 

% 

% This program stores the plot as a metafile with the prefix 'leg' 

% followed by a number starting from 0 on upwards and widing with 
% the extension *.met* You must be familiar with GPP to get hard- 
% copies of the plots. *pllcnl’ is the only variable necessary for 
% the proper indexing of the leg*. met files. Ensure that the 
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% calling program docs not redefine the variable! 

% 

disp([’ ’)’) 

q5=input(’ Store the Plot (Y/N)? 
if q5 == 'V I q5 == ’y\ 
eval([’meta leg\num2str(pltcnt)]) 
pllcnt=pllcnt+l ; 
end 



^*********************** 

% daubdata.m 

% 

% By; CAPT/USMC Kalrabach. M. 

% Professor Lam 

% US Naval Postgraduate School, Monterey CA 

% c>mail: lam@ccc.nps.navy.mil 

% Phone: (408) 646-3044 

% 

% 

function[hcoeff] = daubdata(x) 

% This function returns the proper *h’ coefficients for the specified order 
% Daubechics wavelet. The input value of ’x’ is the specified order, 
if X == 2 

hcoeff = [0.482962913145;0.836516303738;0.224143868042;-0. 129409522551]; 
clscif X = = 3 

hcoeff = [0.332670552950;0.806891509311;0.459877502118;-0. 135011020010; 
-0.085441273882;0.035226291882J; 
clscif X = = 4 

hcoeff = [0.230377813309;0.714846570553;0.630880767930;-0.027983769417; 
-0. 18703481 1719;0.030841381836;0.032883011667;-0. 010597401785]; 
ebeif x = = 5 

hcoeff = [0.160102397974;0.603829269797;0.724308528438;0. 138428145901; 
•0.242294887066;-0.032244869585;0.077571493840;-0. 006241490213; 
-0.012580751999;0.003335725285]; 
ebeif x = = 6 

hcoeff = [0.1 11540743350;0.494623890398;0.751133908021;0. 315250351709; 
-0.226264693965;-0.129766867567;0.097501605587;0.027522865530; 
-0.031582039318;0.000553842201;0.004777257511;-0.001077301085]; 
ebeif x = = 7 

hcoeff = [0.077852054085;0.396539319482;0.729132090846;0.469782287405; 
-0.143906003929;-0.224036184994;0.071309219267;0.080612609151; 
-0.038029936935;-0.016574541631;0.012550998556;0.000429577973; 
-0.001801640704;0.000353713800]; 
ebeif x = = 8 

hcoeff = [0.054415842243;0.312871590914;0.675630736297;0.585354683654; 
-0.01 5829105256;-0.284015542962;0.000472484574;0. 128747426620; 
-0.017369301002;-0.044088253931;0.013981027917;0.008746094047; 
-0.004870352993;-0.000391740373;0.000675449406;-0.0001 17476784]; 
ebeif x = = 9 

hcoeff = [0.038077947364;0.243834674613;0.604823123690;0.657288078051; 
0. 133 197385825;-0.293273783279;-0.096840783223;0. 148540749338; 
0.030725681479;-0.067632829061;0.000250947115;0.022361662124; 
-0.004723204758;-0.004281503682;0.001847646883;0.000230385764; 
-0.000251963189;0.000039347320]; 
ebeif x = = 10 

hcoeff = [0.026670057901 ;0.188176800078;0.527201188932;0. 688459039454; 
0.281 172343661 ;-0.249846424327;-0.195946274377;0. 127369340336; 
0.093057364604;-0.071394147166;-0.029457536822;0.033212674059; 
0.003606553567;-0.010733175483;0.001395351747;0.001992405295; 
-0.000685856695;-0.000116466855;0.000093588670;-0.000013264203]; 
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else 

hcoeff = 0; 
break 
end 

hcoeff = hcoeff /sqrt( 2); % normalize the ’h’ vector 

return 



tf\ /bin/csh -f 
ff 

K metal4 C-Shell Routine 
ff 

# metal4 is a command executable file for the generation of 

# the leg*. met metafiles into postscript graphic files for the 

# SPARC printer number 14 on the fourth floor of Spanagel 431. 

# This routine is system dependent and is guaranteed to work 

# on the UNIX based system of the ECE Department at the Naval 

# Postgraduate School. 

Is -1 leg*. met >tmmmpl 
echo ’#! A)in/csh -f " > tmmmp2 
awk -F. ’{printf "gpp * $1 ' -dps -ol \n* ;\ 
printf "lprl4 * Si ’.ps Vn’;\ 
printf *rm " Si *.ps Vn"}’ tmmmpl > >tmmmp2 
chmod 0700 tmmmp2 
tmmmp2 

rm tmmmpl tmmmp2 
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